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Abstract 

We analyze the signal processing required for the optimal detection of a 
stochastic background of gravitational radiation using laser interferometric 
detectors. Starting with basic assumptions about the statistical properties 
of a stochastic gravity-wave background, we derive expressions for the op- 
timal filter function and signal-to-noise ratio for the cross-correlation of the 
outputs of two gravity- wave detectors. Sensitivity levels required for detec- 
tion are then calculated. Issues related to: (i) calculating the signal-to-noise 
ratio for arbitrarily large stochastic backgrounds, (ii) performing the data 
analysis in the presence of nonstationary detector noise, (iii) combining data 
from multiple detector pairs to increase the sensitivity of a stochastic back- 
ground search, (iv) correlating the outputs of 4 or more detectors, and (v) 
allowing for the possibility of correlated noise in the outputs of two detectors 
are discussed. We briefly describe a computer simulation that was used to 
"experimentally" verify the theoretical calculations derived in the paper, and 
which mimics the generation and detection of a simulated stochastic gravity- 
wave signal in the presence of simulated detector noise. Numerous graphs 
and tables of numerical data for the five major interferometers (LIGO-WA, 
LIGO-LA, VIRGO, GEO-600, and TAMA-300) are also given. This informa- 
tion consists of graphs of the noise power spectra, overlap reduction functions, 
and optimal filter functions; also included are tables of the signal-to-noise ra- 
tios and sensitivity levels for cross-correlation measurements between different 
detector pairs. The treatment given in this paper should be accessible to both 
theorists involved in data analysis and experimentalists involved in detector 
design and data acquisition. 
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I. INTRODUCTION 



The design and construction of a number of new and more sensitive detectors of gravi- 
tational radiation is currently underway. These include the two LIGO detectors being built 
in Hanford, WA and Livingston, LA by a joint Caltech/MIT collaboration the VIRGO 
detector being built near Pisa, Italy by an Italian/French collaboration 0, the GEO-600 
detector being built in Hanover, Germany by an Anglo/German collaboration ||, and the 
TAMA-300 detector being built near Tokyo, Japan M. There are also several resonant bar 
detectors currently in operation, and several more refined bar and interferometric detectors 
presently in the planning and proposal stages. 

The operation of these detectors will have a major impact on the field of gravitational 
physics. For the first time, there will be a significant amount of experimental data to be 
analyzed, and the "ivory tower" relativists will be forced to interact with a broad range of 
experimenters and data analysts to extract the interesting physics from the data stream. 
"Known" sources such as coalescing neutron star (or black hole) binaries, pulsars, super- 
novae, and other periodic and transient (or burst) sources should all be observable with 
gravity-wave detectors. One might also be able to detect a faint stochastic background 
of gravitational radiation, produced very shortly after the big bang. These detections may 
happen soon after the detectors go "on-line" or they may require a decade of further work to 
increase the sensitivity of the instruments. But it is fairly safe to say that eventually, when 
their sensitivity passes some threshold value, the gravity-wave detectors will find sources. 
Even more exciting is the prospect that the detectors will discover new sources of gravi- 
tational radiation — sources which are different from those mentioned above, and which we 
had not expected to find. It promises to be an exciting time. 

The subject of this paper is a stochastic (i.e., random) background of gravitational 
radiation, first studied in detail by Michelson ||, Christensen ||, and Flanagan |7[]. Roughly 
speaking, it is the type of gravitational radiation produced by an extremely large number 
of weak, independent, and unresolved gravity- wave sources. The radiation is stochastic in 
the sense that it can be characterized only statistically. As mentioned above, a stochastic 
background of gravitational radiation might be the result of processes that took place very 
shortly after the big bang. But since we know very little about the state of the universe at 
that time, it is impossible to say with any certainty. A stochastic background of gravitational 
radiation might also arise from more recent processes (e.g., radiation from many unresolved 
binary star systems), and this more recent contribution might overwhelm the parts of the 
background that contain information about the state of the early universe. In any case, the 
properties of the radiation will be very dependent upon the source. For example, one would 
expect a stochastic background of cosmological origin to be highly isotropic, whereas that 
produced by white dwarf binaries in our own galaxy would be highly anisotropic. We will 
just have to wait and see what the detectors reveal before we can decide between these two 
possibilities. 

This paper will focus on issues related to the detection of a stochastic background of 
gravitational radiation. (We will not talk much about possible sources.) We give a complete 
and comprehensive treatment of the problem of detecting a stochastic background, which 
should be accessible to both theorists involved in the data analysis and experimentalists 
involved in detector design and data acquisition. 
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The outline of the paper is as follows: 

In Sec. U we begin by describing the properties of a stochastic background of grav- 
itational radiation — its spectrum, statistical assumptions, and current observational con- 
straints. 

In Sec. PTI] , we describe how one can correlate the outputs of two gravity-wave detectors 
to detect (or put an upper limit on) a stochastic gravity- wave signal. Section [El 1B| includes 
a detailed derivation of the overlap reduction function that covers the case where the two 
arms of a detector are not perpendicular (e.g., GEO-600) and corrects a typographical error 
that appears in the literature. Section [111 U| includes a rigorous derivation of the optimal 
signal processing strategy. Most of the material in Sees. [TT] and [TXT| has already appeared in 
the literature. Interested readers should see Ref. f| for more details, if desired. 

In Sec. [TV], we ask the questions: (i) "How do we decide, from the experimental data, if 
we've detected a stochastic gravity-wave signal?" (ii) "Assuming that a stochastic gravity- 
wave signal is present, how do we estimate its strength?" (iii) "Assuming that a stochastic 
gravity-wave signal is present, what is the minimum value of required to detect it 95% 
of the time?" This leads to a discussion of signal detection, parameter estimation, and 
sensitivity levels for stochastic background searches, adopting a frequentist point of view. 
The calculation of £7q 5%,5% in Sec. |1V D| corrects an error that has appeared in the literature. 

In Sec. |V[ we return to the problem of detection by addressing a series of subtle issues 
initially ignored in Sec. |T|. These include: (i) calculating the the signal-to-noise ratio for 
arbitrarily large stochastic backgrounds, (ii) performing the data analysis in the presence of 
nonstationary detector noise, (iii) combining data from multiple detector pairs to increase 
the sensitivity of a stochastic background search, (iv) correlating the outputs of 4 or more 
detectors, and (v) allowing for the possibility of correlated noise in the outputs of two 
detectors. The material presented in these sections extends the initial treatment of these 
issues given, for example, in Refs. 

Section |V1| consists of a series of graphs and tables of numerical data for the five major 
interferometers (LIGO-WA, LIGO-LA, VIRGO, GEO-600, TAMA-300). The noise power 
spectra, overlap reduction functions, optimal filter functions, signal-to-noise ratios, and sen- 
sitivity levels for cross-correlation measurements between different detector pairs (not just 
LIGO) are given. This information allows us to determine the optimal way of combining 
data from multiple detector pairs when searching for a stochastic background of gravita- 
tional radiation. Section [VI A| also includes a graph of the "enhanced" LIGO detector noise 
curves, which track the projected performance of the LIGO detector design over the next 
decade. 

In Sec. |V1J[ , we describe a computer simulation that mimics the generation and detection 
of a simulated stochastic gravity-wave signal in the presence of simulated detector noise. The 
simulation was used to verify some of the theoretical calculations derived in the previous 
sections. 

Section |VIII| concludes the paper with a brief summary and lists some topics for future 
work. 

Note: Throughout the paper, we use c to denote the speed of light and G to denote 
Newton's gravitational constant (c = 2.998 x 10 10 cm/sec and G = 6.673 x 10~ 8 cm 3 /gm • 
sec 2 ). 
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II. THE STOCHASTIC BACKGROUND - PROPERTIES 



A. Spectrum 

A stochastic background of gravitational radiation is a random gravity-wave signal pro- 
duced by a large number of weak, independent, and unresolved gravity- wave sources. In 
many ways it is analogous to the Cosmic Microwave Background Radiation (CMBR) ||, 
which is a stochastic background of electromagnetic radiation. As with the CMBR, it is 
useful to characterize the spectral properties of the gravitational background by specifying 
how the energy is distributed in frequency. Explicitly, one introduces the dimensionless 
quantity 



^critical am J 



where <ip gw is the energy density of the gravitational radiation contained in the frequency 
range / to / + df, and p C riticai is the critical energy density required (today) to close the 
universe: 

^critical = ^ « 1-6 x 10- 8 h\ m ^ . (2.2) 
87rG cur 



Hq is the Hubble expansion rate (today): 



kin I ! 

H = h m ■ 100 — - = 3.2 x 10- 18 h wo — = 1.1 x 10- 28 c ■ h wo — , (2.3) 

sec • Mpc sec cm 

and hioo is a dimensionless factor, included to account for the different values of Hq that 
are quoted in the literature.^ It is this dimensionless function of frequency, f2 gw (/), that 
we will use to describe the spectrum of a stochastic background of gravitational radiation. 
It follows directly from the above definitions that Q gw (f) h\ m is independent of the actual 
Hubble expansion rate. For this reason, we will often focus attention on this quantity, rather 
than on fl^^f) alone. 

Two remarks are in order: 

(i) There appears to be be some confusion about fl gw (f) in the literature. Some authors 
assume that f2 gw (/) is constant — i.e., independent of frequency. Although this is true for 
some cosmological models, it is not true for all of them. The important point is that any 
spectrum of gravitational radiation can be described by an appropriate f2 gw (/). With the 
correct dependence on frequency, f2 gw (/) can describe a flat spectrum, a blackbody spectrum, 
or any other distribution of energy with frequency. 

(ii) Qgw(f) is the ratio of the stochastic gravity- wave energy density contained in a 
bandwidth Af = f to the total energy density required to close the universe. For the 
CMBR, one can define an analogous quantity: 



hiQQ almost certainly lies within the range 1/2 < h\QQ < 1. 
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Pcritical "111/ 

Since the 2.73 K blackbody spectrum has a peak value of fi e m(/) ^100 ~ 10~ 5 a ^ / = 10 12 
Hz, the CMBR contains (in the vicinity of 10 12 Hz) approximately 10~ 5 of the total energy 
density required to close the universe. A similar interpretation applies to f2 gw (/). 



B. Statistical assumptions 

The spectrum f2 gw (/) completely specifies the stochastic background of gravitational ra- 
diation provided we make enough additional assumptions. We will assume that the stochastic 
background is: (i) isotropic, (ii) unpolarized, (iii) stationary, and (iv) Gaussian. Since these 
properties might not hold in general, it is worthwhile to consider each one of them in turn. 

(i) Since it is now well established that the CMBR is highly isotropic [Q, it is not unrea- 
sonable to assume that a stochastic background of gravitational radiation is also isotropic. 
But this assumption might not be true. For example, as mentioned in Sec. |, if the domi- 
nant source of the stochastic gravity-wave background is a large number of unresolved white 
dwarf binary star systems within our own galaxy, then the stochastic background will have 
a distinctly anisotropic distribution, which forms a "band in the sky" distributed roughly 
in the same way as the Milky Way galaxy. It is also possible for a stochastic gravity-wave 
background of cosmological origin to be anisotropic, although one would then have to ex- 
plain why the CMBR is isotropic but the gravity- wave background is not. In either case, 
such anisotropies can be searched for in the data stream. (See Ref. []TD for details.) 



(ii) The second assumption is that the stochastic gravity-wave background is unpolarized. 
This means that the gravitational radiation incident on a detector has statistically equivalent 
"plus" and "cross" polarization components. We see no strong reason why this should not 
be the case. 

(iii) The assumption that the stochastic background is stationary (i.e., that all statistical 
quantities depend only upon the difference between times, and not on the choice of time 
origin) is almost certainly justified. This is because the age of the universe is roughly 20 
orders of magnitude larger than the characteristic period of the waves that LIGO, VIRGO, 
etc. can detect, and 9 orders of magnitude larger than the longest realistic observation times. 
It seems very unlikely that a stochastic background of gravitational radiation would have 
statistical properties that vary over either of these time-scales. But unlike the stochastic 
gravity-wave background, the noise intrinsic to the detectors will change over the course of 
the observation times. This poses a problem for the data analysis, which we initially ignore 



in Sec. |T|. We return to this problem in Sec. |VB| where we discuss nonstationary detector 
noise. 

(iv) The final assumption is that the stochastic gravity-wave background is a Gaussian 
random process. This means that the joint probability density function of the gravitational 
strains hifa), hj(tj), ■ ■ ■ in detectors i, j, ■ ■ ■ is a multivariate Gaussian (i.e., normal) dis- 
tribution. In this case, the mean values (Kit)) and the second-order moments {hi(ti)hj(tj)) 
completely specify the statistical properties of the signal. For many early-universe processes, 
or even for more recent sources of a gravity-wave background, this is a reasonable assump- 
tion. It can be justified by the central limit theorem, which says that the sum of a large 



5 



number of statistically independent random variables is a Gaussian random variable, inde- 
pendent of the probability distributions of the original variables. This will be the case for 
the stochastic background if it is the sum of gravity-wave signals produced by a large num- 
ber of independent gravity-wave sources. This assumption will not be true, however, if the 
stochastic background is the sum of the radiation produced, e.g., by only a few unresolved 
binary star systems radiating in a given frequency interval at any instant of time. (See, e.g., 
Ref. [0.) 

The above four properties form the basis for the statistical analysis that we will give in 
the following sections. We will assume that they hold throughout, unless we explicitly state 
otherwise. 



C. Expectation value 

Using the definition of the spectrum Q gw (/) and the statistical assumptions described in 
the previous subsection, we can derive a useful result for the expectation value of the Fourier 
amplitudes of a stochastic background of gravitational radiation. This result will be needed 
in Sec. |T| when we discuss signal detection and optimal filtering. 

The starting point of the derivation is a plane wave expansion for the gravitational metric 
perturbations in a transverse, traceless gauge: 

h ab (t, x) = £ /°° df [ dCl h A (f, fl) e °^M-*/«) 4(0) . (2.5) 

Here O is a unit vector specifying a direction on the two-sphere, with wave vector k := 
2-nffl/c. Also, e^ b (fl) are the spin-two polarization tensors for the "plus" and "cross" po- 
larization states A = +, x. Explicitly, 

e+ fe (0) = rh a m b - h a h b , (2.6) 

e* b (Vt) = rh a h b + h a rh b , (2.7) 

where 

Cl = cos sin 9 x + sin sin 9 y + cos 9 z , (2-8) 

rh = sin <f) x — cos cf) y , (2-9) 

h = cos cos 9 x + sin cos 9 y — sin 9 z , (2-10) 

and (9, 0) are the standard polar and azimuthal angles on the two-sphere. The Fourier 
amplitudes ^(Z, O) are arbitrary complex functions that satisfy Ha{— f, ^) = h* A (f,fl), 
where * denotes complex conjugation. This last relation follows as a consequence of the 
reality of h ab (t, x). 

The assumptions that the stochastic background is isotropic, unpolarized, and stationary 
imply that the expectation value (i.e., ensemble average) of the Fourier amplitudes ^(/, O) 
satisfies: 

(h A (f,n)h A/ (f,Cl')) = 5 2 (h,Cl')5 AA> 5(f - f) H(f) , (2.11) 
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where 5 2 (Cl,Cl') := 5(<f) — <j)')5(cos6 — cos#') is the covariant Dirac delta function on the 
two-sphere, and H(f) is a real, non-negative function, satisfying H(f) = H(—f)Q If we 
further assume that the stochastic background has zero mean, then 

(h A (f,Q)) = 0. (2.12) 

Finally, since we are assuming that the stochastic background is Gaussian, the expectation 
values (|2.11| ) and (|2.12|) completely specify its statistical properties. 

H(f) is related to the spectrum Q gw (/) of the stochastic gravity-wave background. This 
follows from the expression 

Pgw = ^ (h ab (t,x)h ab (t,x)) (2.13) 



for the energy density in gravitational waves (see, e.g., p. 955 of Ref. |12[). By differentiating 
the plane wave expansion ( |2.5|) with respect to t, forming the contraction in Eq. (|2.13|) , and 
calculating the expectation value using ( |2.12| ), we find 

ft. = nrf*' ,jr </> (==f*^)- (2 - 14) 

Using Eqs. Q and for Q gw {f), then yields 



3/7 2 

H(f) = 3^ I/I" 3 OgwCI/l) • (2-15) 



Thus, 



(h* A (f,{i)h AI (f',&)) = |S ^^ftoMz-zo i/r 3 fi S w(i/i) , (2.i6) 

which is the desired result. 



D. Observational constraints 

At present, there are three observational constraints on the stochastic gravity-wave spec- 
trum flgw(f). These constraints are quite weak in the frequency range of interest for 
ground-based interferometers (1 Hz < / < 10 3 Hz) and for proposed space-based detec- 
tors (10 -4 Hz < / < 10 _1 Hz). There are tighter constraints on the spectrum in two 
frequency ranges, and one "wideband" but very weak constraint. In this paper, we sim- 
ply state the constraints. For a more complete discussion, see Ref. |§ and the references 
mentioned therein. 



2 If the stochastic background is anisotropic, we should replace H(f) by a function that depends 
on Ct in addition to /. If the stochastic background is polarized, we should replace H(f) by a 
function that depends on the polarization A = +, x as well. 
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(i) The strongest observational constraint on Q gw (f) comes from the high degree of 



isotropy observed in the CMBR. In particular, the one- year [|13yi4l , two-year [15|], and four- 
year |16j data sets from the Cosmic Background Explorer (COBE) satellite place very strong 
restrictions on fl^f) at very low frequencies: 




flgw(/) K 00 < 7 x lO" 11 -H for H <f<30H . (2.17) 



Note that the above constraint does not apply to any gravitational wave, but only to those 
of cosmological origin that were already present at the time of last scattering of the CMBR. 
Also, since H = 3.2 x 10~ 18 h WQ Hz, this limit applies only over a narrow band of frequencies 
(10~ 18 Hz < / < 10~ 16 Hz), which is far below any frequency band accessible to investigation 
by either earth-based or space-based interferometers. Thus, although this constraint is 
severe, it is not directly relevant for any of the present-day gravity-wave experiments. 

(ii) The second observational constraint comes from almost a decade of monitoring the 



radio pulses arriving from a number of stable millisecond pulsars [17]. These pulsars are 
remarkably stable clocks, and the regularity of their pulses places tight constraints on f2 gw (/) 
at frequencies on the order of the inverse of the observation time of the pulsars (~ Hz): 

figwCf = KT 8 Hz) h 2 100 < 1(T 8 . (2.18) 

Like the constraint on the stochastic gravity-wave background from the isotropy of the 
CMBR, the millisecond pulsar timing constraint is irrelevant for current gravity-wave ex- 
periments. The frequency / = 10~ 8 Hz is 10 orders of magnitude smaller than the band of 
frequencies accessible to LIGO, VIRGO, etc., and 4 orders of magnitude smaller than that 
for proposed space-based detectors. 

(hi) The third and final observational constraint on Q gw (/) comes from the standard 
model of big-bang nucleosynthesis || . This model provides remarkably accurate fits to the 
observed abundances of the light elements in the universe, tightly constraining a number 
of key cosmological parameters. One of the parameters constrained in this way is the ex- 
pansion rate of the universe at the time of nucleosynthesis. This places a constraint on the 
energy density of the universe at that time, which in turn constrains the energy density in 
a cosmological background of gravitational radiation: 



ciln/^ gw (/)/ ioo <10- 5 . (2.19) 

/>10- 8 Hz 

Although this bound constrains the spectrum of gravitational radiation fl^i^f) over a broad 
range of frequencies, it is not very restrictive. 



III. THE STOCHASTIC BACKGROUND - DETECTION 

In this section, we begin our detailed discussion of the detection of a stochastic back- 
ground of gravitational radiation. We explain how one can correlate the outputs of two 
gravity-wave detectors to detect (or put an upper limit on) a stochastic background signal. 
In Sec. |111 B|, we give a detailed derivation of the overlap reduction function that covers the 
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case where the two arms of a detector are not perpendicular (e.g., GEO-600). In Sec. |III C . 
we give a rigorous derivation of the optimal signal processing strategy. The statistical as- 
sumptions that we will use for the stochastic gravity-wave background are those described 
in Sec. [II B| . In addition, we will assume that the noise intrinsic to the detectors are: (i) 
stationary, (ii) Gaussian, (iii) statistically independent of one another and of the stochastic 
gravity-wave background, and (iv) much larger in magnitude than the stochastic gravity- 
wave background. The modifications that are necessary when one relaxes most of these 
assumptions will be discussed in Sec. |V|. 



A. Coincident and coaligned detectors 

To begin, let us consider the simplest possible case. Let us suppose that we have two 
coincident and coaligned gravity-wave detectors with outputs 

si(t) := h(t) + m(t) , (3.1) 
s 2 {t) := h 2 {t) +n 2 {t) . (3.2) 

Here hi(t) and h 2 (t) denote the gravitational strains in the two detectors due to the stochastic 
background, and ni{t) and n 2 (t) denote the noise intrinsic to the first and second detector, 
respectively.^ Since we are assuming that the two detectors are coincident and coaligned 
(i.e., have identical locations and arm orientations), the gravitational strains are identical: 

h{t) := h^t) = h 2 {t) . (3.3) 

But the noise ni(t) and n 2 (t) are not equal to one another. As mentioned above, we will 
assume that they are stationary, Gaussian, statistically independent of one another and of 
the gravitational strain, and much larger in magnitude than the gravitational strains.^ 

Given the detector outputs Si(t) and s 2 (t), we can form a product "signal" 5* by multi- 
plying them together and integrating over time: 

,T/2 

S := / dt Sl (t)s 2 (t) . (3.4) 

J-T/2 

This quantity is proportional to the (zero-lag) cross-correlation of si(t) and s 2 (t) for an 
observation time T. Since si(t) and s 2 (t) are random variables, so too is S. It has a mean 
value 

A* := (S) , (3.5) 



3 We will assume throughout that the detector outputs are not whitened. 

4 The assumption that the noise intrinsic to the detectors are statistically independent of one 
another is unrealistic for the case of coincident and coaligned detectors. But it is a reasonable 
assumption for widely-separated detector sites. (See Sec. VE for more details.) 
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and variance 

a 2 := (S 2 ) - (S) 2 , (3.6) 

which are related to the variances of rii(t), n 2 (t), and The goal is to calculate \i and 

a, and then to construct the signal-to-noise ratio 

SNR := ^ . (3.7) 
a 

As we shall see in Sec. |V|, the value of the signal-to-noise ratio enters the decision rule for 
the detection of a stochastic gravity-wave signal. 
Let us start with the mean value \i. By definition, 

rT/2 

H ■= (S) = / dt (s 1 (t)s 2 (t)) (3.8) 

J-T/2 

r T/2 

dt (h 2 (t) + h(t)n 2 (t) + nx(t)h(t) + n x {t)n 2 (t)) (3.9) 

-T/2 
r T/2 

dt (h 2 (t)) (3.10) 

-T/2 

= T (h 2 (t)) =:To* h , (3.11) 

where o\ denotes the (time- independent) variance of Note that we used the statistical 

independence of ni(t), n 2 (t), and h(t) to obtain the third line, and the stationarity of h(t) 
to obtain the last. 

To express the variance a\ := (h 2 (t)) in terms of the frequency spectrum fl^^f), we will 
make use of the plane wave expansion ( |2.5| ) and the expectation value ( |2.16| ). Since 

h{t) := h ab (t,x ) l - (X a X b - Y a Y b ) (3.12) 

(where Xq is the common position vector of the central station of the two coincident and 
coaligned detectors, and X a and Y a are unit vectors pointing in the directions of the detector 
arms),[] it follows that 



5 The mean values of rii(t), n 2 {t)^ and h(t) are equal to zero, either by assumption or by definition. 

6 The dimensions of a\ and a 2 are different: a\ has dimensions of strain 2 , while a 2 has dimensions 
of strain 4 • sec 2 . (See, e.g., Eq. (^20|) .) 

7 xq and X a , Y a are actually functions of time due to the Earth's motion with respect to the 
cosmological rest frame. They can be treated as constants, however, since: (i) the velocity of the 
Earth with respect to the cosmological rest frame is small compared to the speed of light, and (ii) 
the distance that the central stations and arms move during the correlation time between the two 
detectors is small compared to the arm length. (The correlation time equals zero for coincident and 
coaligned detectors; it equals the light travel time between the two detectors when the detectors 
are spatially separated.) See Ref. [1C] for more details. 
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" df \f\~ 3 n^d/l) , (3.16) 



oo 



r f ~ rinfty poo 

^ = EE/ dn dn' df df (h* A (f,n)h A ,(f,n')) 

A A' J— oo J— oo 

x e -itorf(t-fr3 /c) e i2nf(t-Q'-x /c) pAfypA' (ft} ? (3^3) 

where 

:= e^ 5 (fi) i (x a X 6 - F a r 6 ) (3.14) 

is the response of either detector to a zero frequency, unit amplitude, A = +, x polarized 
gravitational wave. Using ( |2.16| ) for the expectation value {h A (f,Cl)h A '(f',Cl')), the above 
expression for a 2 simplifies to 

= # I/!" 3 E / 52 dn F A (Q)F A (Q) (3.15) 

3# p 2 r 00 

~ 20tt 
where we used 

Y.J S2 dVF A (n)F A (£l) = ^ (3.17) 

to obtain the last line. Thus, for coincident and coaligned detectors, the mean value of the 
cross-correlation signal S is 

"iff 2 roo 

» = ^ T df |/|- 3 n^d/l) . (3.18) 

This is the first of our desired results. 

To evaluate the variance a 2 , we will make use of the assumption that the noise intrinsic 
to the detectors are much larger in magnitude than the gravitational strain. Then 

a 2 := (S 2 ) - {S} 2 « (S 2 ) (3.19) 

fT/2 fT/2 

= dt dt' (s 1 (t)s 2 {t)s 1 {t')s 2 (t')) (3.20) 

J-T/2 J-T/2 
T/2 pT/2 

dt / dt' (ni(*)n2(t)ni(On 2 (0) ( 3 - 21 ) 

-T/2 J-T/2 
T/2 pT/2 

dt / rft' (n 1 (t)n 1 (t'))(n 2 (t)n 2 (t')) , (3.22) 

-T/2 J-T/2 

where we used the statistical independence of n\{t) and n 2 (t) to obtain the last line. By 
definition^] 



8 Equation ( 3.23| ) can also be written in the frequency domain: 

(n*(f)Mf)) = ^(f-f)P t (\f\). 
See the discussion surrounding Eq. ( |3.64| ) for more details. 
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A J — oo 



df ^ i{t - l,) Pi(\f\) 



(3.23) 



where Pi(\f\) is the (one-sided) noise power spectrum of the ith detector (i = 1,2). Pi(\f\) 
is a real, non-negative function, defined with a factor of 1/2 to agree with the standard 
(one-sided) definition used by instrument-builders. It satisfies^ 



<"?(«)> 



df P(f) 



(3.24) 



so the total noise power is the integral of Pi(f) over all positive frequencies / from to oo, 
not from — oo to oo. (Hence the reason for the name one-sided.) Graphs of the predicted 
noise power spectra for the initial and advanced LIGO detectors are shown in Fig. [TJ. Graphs 
of the predicted noise power spectra for the other major interferometers (i.e., VIRGO, GEO- 
600, and TAMA-300) and for the "enhanced" LIGO detectors are shown in Figs. |TT| - [H)| in 



Sec. VIA 



LIGO Noise Power Spectra 




Advanced LIGO 

1000 10000 



f (Hz) 



FIG. 1. A log-log plot of the predicted noise power spectra for the initial and advanced LIGO 
detectors. The data for these noise power spectra were taken from the published design goals fij. 



Inserting (|3~23p into (ggg) yields 



cr 2 ps 



■T/2 ,T/2 

dt / dt! 

4 J -T/2 J -T/2 ./-x 



df / df 

J —oo 

x e^/C-O e -^f'it-t') Pl (\f\)p 2 (\f\) s 



(3.25) 



where we used the reality of n^if) and P2(\f\) to produce the minus sign in the power of the 
second exponential. If we integrate this expression over t and t', we find 



9 Unlike a\ and a 2 , a\ and a 2 , have the same dimensions (strain 2 ). 
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1 poo poo 

o 2 ^-A df df'SUf-f')Pi(\f\)P2(\f'\), (3-26) 

*± J —oo J —oo 

where 

M/) := r dt e*** = Sin(?r / T) (3.27) 

J-r/2 7T / 

is a finite-time approximation to the Dirac delta function 5(/). In the limit T — > oo, 
&r(/) reduces to 5(f), but for a finite observation time T, one has 5t(0) = T. Since in 
practice the observation time T will be large enough so that 5x{f — /') is sharply peaked 
over a region in / — /' whose size ~ 1/T is very small compared to the scale very on which 
the functions Pi(|/|) and Ad/I) are varying]^ we can replace one of the finite-time delta 
functions 8r{f — f) by an ordinary Dirac delta function, and evaluate the other at / = /'. 
Doing this yields 

° 2 ~-A df P^fDP.dfl) , (3.28) 

4 J —oo 

which is the second of our desired results. 

Using ( |3.18| ) and ( |3.28j ), we can now form the signal-to- noise ratioQ 



M^.Myf !-^f\f\'^M) (329) 

a Wt,- r r rf/ ^(1/1)^(1/1)1 



The multiplicative factor of yT means that we can always exceed any prescribed value 
of the signal-to-noise ratio by correlating the outputs of two gravity-wave detectors for a 
long enough period of time.0 We will have more to say about signal detection, parameter 
estimation, and sensitivity levels for stochastic background searches in Sec. [TV] . 



B. Overlap reduction function 

To provide a rigorous treatment of the signal analysis for a stochastic background of 
gravitational radiation, we must take into account the fact that the two gravity-wave de- 
tectors will not necessarily be either coincident or coaligned. There will be a reduction in 



10 Typically, an observation time T will be on the order of months (i.e., 10 7 seconds), while the 
noise power spectra Pi{\f\) vary on a scale of greater than a few Hz. 



11 Warning!! The signal processing strategy described above, leading to Eq. ( p. 29 ), is not optimal 
even for the case of coincident and coaligned detectors. The optimal signal processing strategy, 



which is described in Sec. QIC , leads to the signal-to- noise ratio given by Eq. (|77|). Setting 
■y(f) = 1 in Eq. ( p. 75 ) yields an expression for the optimally-filtered signal-to-noise ratio for the 
case of coincident and coaligned detectors. 



12 This assumes that there is no systematic source of correlated detector noise. In Sec. VE, we 
discuss the limits that correlated detector noise impose. 
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sensitivity due to: (i) the separation time delay between the two detectors, and (ii) the 
non-parallel alignment of the detector arms. These two effects imply that h\(t) and /i 2 (^) 
are no longer equal; the overlap between the gravitational strains in the two detectors is 
only partial. Statistically, these effects are most apparent in the frequency domain. 

The overlap reduction function 7(/), first calculated in closed-form by Flanagan 
quantifies these two effects. This is a dimensionless function of frequency /, which is deter- 
mined by the relative positions and orientations of a pair of detectors. Explicitly, 




(3.30) 



where Cl is a unit vector specifying a direction on the two-sphere, Ax := x\ — £2 is the 
separation vector between the central stations of the two detector sites, and 

:= e£(fi) df := ei(fl) ~ - Y?Y*) (3.31) 

is the response of the ith detector (i = 1,2) to the A = +, x polarization. (See also 
Eq. ( p.!4| ).) The symmetric, tracefree tensor df 3 specifies the orientation of the two arms 
of the ith detector. The overlap reduction function 7(/) equals unity for coincident and 
coaligned detectors. It decreases below unity when the detectors are shifted apart (so there 
is a phase shift between the signals in the two detectors), or rotated out of coalignment (so 
the detectors are sensitive to different polarizations). In Sec. [Ill C| , we will see that j(f) 
arises naturally when evaluating the expectation value of the product of the gravitational 
strains at two different detectors when they are driven by an isotropic and unpolarized 
stochastic background of gravitational radiation. 

To get a better feeling for the meaning of 7(/), let us look at each term in Eq. (|3.30| ) 
separately: (i) The overall normalization factor 5/87T is chosen so that for a pair of coincident 
and coaligned detectors 7(/) = 1 for all frequencies /. (ii) The sum over polarizations A is 
appropriate for an unpolarized stochastic background, (iii) The integral over the two-sphere 
is an isotropic average over all directions Cl of the incoming radiation, (iv) The exponential 
phase factor is the phase shift arising from the time delay between the two detectors for 
radiation arriving along the direction Cl. In the limit / — ► 0, this phase shift also goes to 
zero, and the two detectors become effectively coincident, (v) The quantity J2a FfityF^ity 
is the sum of products of the responses of the two detectors to the + and x polarization 
waves. For coaligned detectors, Ff(Q) = F 2 A (f2) and the integral of this quantity over the 
two-sphere equals the inverse of the overall normalization factor. (See Eq. (|3.17f) .) 

Figure |2] shows a graph of the overlap reduction function 7(/) for the Hanford, WA 
and Livingston, LA LIGO detector pair.0 Note that the overlap reduction function for the 
LIGO detector pair is negative as / — > 0. This is because the arm orientations of the two 
LIGO detectors are not parallel to one another, but are rotated by 90°. If, for example, the 
Livingston, LA detector arms were rotated by 90° in the clockwise direction, only the overall 



b igures in Sec. |VIB| show graphs of the overlap reduction functions for different detector 

pairs. 
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sign of 7(/) would change. Note also that the magnitude of 7(0) is not unity, because the 
planes of the Hanford, WA and Livingston, LA detectors are not identical.^ Thus, the arms 
of the two detectors are not exactly parallel, and (7(0) | = 0.89, which is less than 1. 



0.2 

0.0 
-0.2 
^ -0.4 
-0.6 
-0.8 
-1.0 

0.0 100.0 200.0 1 10 100 1000 10000 

f (Hz) f (Hz) 

FIG. 2. The overlap reduction function 7(/) for the Hanford, WA and Livingston, LA LIGO 
detector pair. (The horizontal axis of the left-hand graph is linear, while that of the right-hand 
graph is logio-) The overlap reduction function has its first zero at 64 Hz, as explained in the text. 
It falls off rapidly at higher frequencies. 

From Fig. [|, one also sees that the overlap reduction function for the two LIGO detectors 
has its first zero at 64 Hz. This can be explained by the fact that a gravitational plane wave 
passing by the earth excites a pair of detectors in coincidence when the positive (or negative) 
amplitude part of the wave is passing by both detectors at the same time; it excites the two 
detectors in anti-coincidence when the positive (or negative) amplitude part of the wave 
is passing by one detector, and the negative (or positive) amplitude part of the wave is 
passing by the other detector. (See Fig. [|) Provided that the wavelength of the incident 
gravitational wave is larger than twice the distance between the two detectors, the detectors 
will be driven in coincidence (on the average). For the case of the LIGO detector pair, this 
means that the Hanford, WA and Livingston, LA detectors will be driven in coincidence 
(on the average) by an isotropic and unpolarized stochastic background of gravitational 
radiation having a frequency of less than / = c/(2d) = 50 Hz. The actual frequency of the 
zero (/ = 64 Hz) is slightly larger than this, since 7(/) is a sum of three spherical Bessel 
functions, which does not vanish at exactly 50 Hz. 

In Appendix B of Ref. J7], Flanagan outlines a derivation of a closed-form expression 
for the overlap reduction function 7(/). The resulting expression applies to any pair of 
gravity- wave detectors, including interferometers with nonperpendicular arms and/or arbi- 
trary orientations. This is a useful result, because, e.g., the arms of the GEO-600 detector 



The two LIGO detectors are separated by an angle of 27.2° as seen from the center of the earth. 
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Min Zero Max 




FIG. 3. The surface of the earth (15° < latitude < 75°, -130° < longitude < 20°) including 
the LIGO detectors in Hanford, WA (LI) and Livingston, LA (L2), the VIRGO detector (V) in 
Pisa, Italy, and the GEO-600 (G) detector in Hanover, Germany. The perpendicular arms of the 
LIGO detectors are also illustrated (though not to scale). A plane gravitational wave passing by 
the earth is indicated by successive minimum, zero, and maximum of the wave. As this wave 
passes by the pair of LIGO detectors, it excites the two in coincidence at the moment shown, since 
both detectors are driven negative by the wave. During the time when the zero is between LI and 
L2, the two detectors respond in anti-coincidence. Provided that the wavelength of the incident 
gravitational wave is larger than twice the separation (d = 3001 km) between the detectors, the 
two detectors are driven in coincidence more of the time than in anti-coincidence. 



are separated by 94.33°. Below we give a more detailed version of the derivation that appears 
in Ref. M, and correct a typographical error that appears in Eq. (B6) of that paper. 

We take, as our starting point for the derivation, the integral expression (|3.30| ) for 7(/). 
To simplify the notation in what follows, we also define 

Ax := d s and a := — — , (3.32) 

c 

where s is a unit vector that points in the direction connecting the two detectors, and d is 
the distance between the two detectors. In terms of these quantities, we can write 

1 {f) = dfd< 2 d T ahcd {a 1 s) , (3.33) 

where 

W", S) :=^£ / dCl e iatl - s e£(n)e£(n) . (3.34) 
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r a fecd(«, s) is a tensor which is symmetric under the interchanges a <-> 6, c <-> d, afe <-> erf. It 
is also tracefree with respect to the a& and cc? index pairs. 

To evaluate T abcd (a,s), we begin by writing down the most general tensor constructed 
from 5 ab and s a that has the above-mentioned symmetry properties: 



r a fe C d(a, s) = A(a) 5 ab 5 cd + B(a) {5 ac 5 bd + 5 bc 5 ad ) +C(a) {5 ab s c s d + 5 cd s a s b ) 
+D(a) (5 ac s b s d + 5 ad s b s c + 5 bc s a s d + 5 bd s a s c ) + E(a) s a s b s c s d 



(3.35) 



We then contract flOgD with 5 ab 5 cd , (5 ac 5 bd + 5 bc 5 ad ), 
of equations for the functions A, B , • ■ ■ , E: 



s a s s c s to obtain a linear system 



where 
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p(a) 
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[a, 


3) 


5 ab 5 cd ^ 


g(a) 


'■= Fabcd 


[a, 


s) 


(5 ac 5 bd + 6 bc 6 ad ) , 


r(a) 


'■— Fabcd 


[a, 


s) 


(S ab s c s d + 5 cd s a s b ) , 


s(a) 


'■= F a bcd 


[a, 


s) 


(5 ac s b s d + 5 ad s b s c + 5 bc s a s d + 5 bd s a s c 


t(a) 


'■= Tafecrf 


>, 


s) 


s a s b s c s d . 



(3.36) 



(3.37) 



From Eq. G3.34Q , we see that the functions p,q,---,t are scalar integrals that involve con- 
tractions of the spin- two polarization tensors ef b (Cl). 

To evaluate these integrals, we choose (without loss of generality) a coordinate system 
where the unit vector s coincides with unit vector z. Then 



and 



s = 


cos 6 , 


rh ■ 


p(a 


)=o, 




q(a 


) = 10 


1 dx 


r{a 


) = o, 




s(a 


) = 10 


1 dx 


t(a 




dx 



, n ■ s 



sin 9 , 



20 j (a) 



e iax ^ 



X 



„2\2 



X 



40 

— .71(a) , 
a 

20 

— 32(a) , 
cr 



where jo(a), ji(oc), and .72 (a) are the standard spherical Bessel functions: 



(3.38) 



(3.39) 
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Jo(aJ 
Ji(a) 
32(a) 



sin a 

a 
sin a 



a- 



cos a 
a 



(3.40) 



sm a cos a 
3 — =--3 — =- 



or 



a- 



sma 
a 



Note that p(a) = and r(a) = are immediate consequences of the tracefree property of 

T abcd (a,s). 

The above linear system of equations ( 3.36|) can be inverted for the functions A,B, - ■ ■ ,E. 
The results are: 



'A' 




' -ha 2 


10a 


5" 




B 




ha 2 


-10a 


5 




' Jo 


C 
D 


ha 2 
-ha 2 


-10a 
20a 


-25 
-25 




Ji 

J2 


E 




ha 2 


-50a 


175 





a 



(3.41) 



Finally, to obtain an expression for the overlap reduction function 7(/), we substitute (|3.35|) 



into ( |3.33| ). Since d. 



ab 



1, 2) is tracefree, it follows that 



7 (/) = 25(a) df d 2ab + AD(a) dfd 2a c s b s c + E(a) dfd c 2 d s a s b s c s d . (3.42) 
Substituting the expressions for the functions B, D, E given by (|3.41| ) into ( |3.42|) yields 



7(/) = PM) dfd 2a b + /02(a) dfd 2a c s b s c + p 3 (a) dfdfs a s h s c s d 



where 



' Pi ' 




P2 




.P3_ 





10a 2 -20a 10 
-20a 2 80a -100 
5a 2 -50a 175 



This is the desired result. 

Note: In Eq. (B6) of Ref. 0, the factor n 
in Eq. ( |3.44| ), this factor should equal —10/ a. 





'jo' 




jl 




J2 . 


Jl 


(a) i 



a 



(3.43) 



(3.44) 



C. Optimal filtering 

Using the techniques developed in the previous two subsections, we are now in a position 
to give a rigorous derivation of the optimal signal processing required for the detection of a 
stochastic background of gravitational radiation. We start by writing the cross-correlation 
signal S between the outputs of the two detectors in the following form: 



S :-- 



T/2 
-T/2 



dt 



T/2 
-T/2 



dt' Sl (t)s 2 (t')Q(t,t') , 



(3.45) 



where as before 
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Sl (t) := h^+rnit) , 
s 2 (t) := h 2 {t)+n 2 (t) , 



(3.46) 
(3.47) 



but now Q(t, t') is a filter function, which is not necessarily equal to 5(t — t') as we assumed 
in Sec. |111 A| . Because we are assuming in this section that the statistical properties of the 
stochastic gravity-wave background and noise intrinsic to the detectors are both stationary, 
the best choice of filter function Q(t, t') can depend only upon the time difference At := t— t' . 
The goal is to find the optimal choice of filter function Q(t — t') := Q(t, t') in a rigorous way. 

The optimal choice of filter function Q(t — t') will depend upon the locations and orienta- 
tions of the detectors, as well as on the spectrum of the stochastic gravity-wave background 
and the noise power spectra of the detectors. It falls off rapidly to zero for time delays 
At = t — t' whose magnitude is large compared to the light travel time d/c between the two 
sites.0 (See Fig. ||, which is located at the end of this section.) Since a typical observation 
time T will be ^> d/c, we are justified in changing the limits on one of the integrations in 
( P.45| ) to obtain 

i-T/2 i-oo 

S= dt dt' s l (t)s 2 (t')Q(t - 1') . (3.48) 

J-T/2 J-oo 

This change of limits simplifies the mathematical analysis that follows. 

We can also write Eq. (|3.48| ) in the frequency domain. Using the convention 

/oo 
dt c-^ g(t) (3.49) 
-oo 

for the Fourier transform of g(t), it follows that 

/OO fOQ _ 

df / df'5 T (f-f')sl(f)s 2 (f)Q(f'), (3-50) 
-oo J — oo 

where 5i(/), s 2 (/), and Q(f) are the Fourier transforms of Si(t), s 2 (t), and Q(t — t'), and 
<5r(/ — /') is t ne finite-time approximation to the Dirac delta function 5(f — /') defined by 
( Pn?7| ). Note also that for a real Q(t - t'), Q(-f) = Q*(f). 

The optimal choice of filter function also depends on the quantity that we want to 
maximize. As we shall see in Sec. [TV], in the context of stochastic background searches, it is 
natural to maximize the signal-to-noise ratio 

SNR := - , (3.51) 
a 

where [i and a 2 are the mean value and variance of the cross-correlation signal S defined 
by ( p.5[ ) and (|3~B|). The techniques that we will use to evaluate /i and a 2 are very similar to 
those that we used in Sec. [Ill A| for the case of coincident and coaligned detectors. 

The calculation of the mean value /i is straightforward. Since we are assuming that the 
noise intrinsic to the two detectors are statistically independent of each other and of the 
gravitational strains, it follows immediately from (|3.50|) that 



d/c = 10 2 sec for the LIGO detector pair. 
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fi := (S) = df df 8 T (f - f')(K(f)h2(f))Q(f) • (3-52) 

J —oo J —oo 

To calculate the expectation value (h*(f) hi (/'))> we again make use of the plane wave 
expansion ( ^5) and the expectation value (|2.16|) . Since 

Hf) = E [dti h A (f,Q) e- i2 ^/ c F^Q) , (3.53) 
where % — 1, 2 labels the two detectors, we find 

(^(/)M/')) = EE l/h f dh' (h* A (f,h)h A if'M 

A A' S S 

x j2*f€l.x t /c e -i2Kf>fl>-x 2 /c F f^)Ff{Cl') (3.54) 

<k/-/o i/r 3 fi gw (i/D 



32vr 3 

x 



E ^ 2 d£l e t2nfCl - A5;/c Ff((l)Ff((l) (3.55) 



Mi 
20vr 2 



*(/-/0 l/r 3 ^gw(|/|)7(l/l), (3.56) 



where we used Q2.16D to obtain the second equality and the definition ( |3.30| ) of the overlap 
reduction function to obtain the third. Substituting (|3.56|) into ( |3.52| ) yields 

"iff 2 r oo 

»=^ T j_Jf l/l" 3 fi g w(|/|) 7(1/1) W) • (3.57) 

The factor of T on the RHS arises from evaluating St{0)- 

Before calculating the variance a 2 , it is worthwhile to make a slight digression and study 
in more detail the expectation value ( |3.56|) derived above. It turns out that this equation 
has an important physical implication. In terms of the time domain variables hi(t) and 
hiit'), Eq. ( p.56| ) can be rewritten as 

/oo 
df e^*-*) H a (f) , (3.58) 
-oo 

where 

Hu(J) = ^ l/l" 3 n^d/l) 7(1/1) • (3.59) 

In other words, Hiz(f) is just the Fourier transform of the cross-correlation of the gravita- 
tional strains hi(t) and h2{t') at the two detector sites. Moreover, since the noise intrinsic 
to the two detectors are statistically independent of one another and of the gravitational 
strains, it follows that 

(MtWO) = < s iW s 2(0> • (3.60) 

Thus, Hi 2 (f) is the Fourier transform of the cross-correlation of the outputs of the two 
detectors. But this correlation is something that we can measure (or at least estimate) given 
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enough data.0 This in turn implies that we can measure (or at least estimate) f2 gw (|/|). 
Explicitly, given the measured values of {si(t)s2(t')), we take their Fourier transform (with 
respect t — t'), multiply by |/| 3 , and divide by 7(|/|), to determine fi gw (|/|). This will 
yield a good approximation to the real stochastic gravity-wave spectrum provided that the 
noise intrinsic to the detectors is not too large, or equivalently, if we measure the detector 
outputs for a long enough period of time T. Also, the approximation of £\w(|/|) will be 
best for frequencies l/T < f < c/d (where d/c is the light travel time between the two 
detectors): For / < l/T, the Fourier transform lacks sufficient time domain data to provide 
useful information; for f > c/d, the overlap reduction function quickly approaches zero, so 
the division by 7(/) makes f^ gw (|/|) ill-behaved. 

Let us return now to the calculation of the signal-to- noise ratio SNR := fi/a. To calculate 
the variance a 2 , we assume (as in Sec. pTA|) that the noise intrinsic to the two detectors are 



much larger in magnitude than the stochastic gravity-wave background. Then 

a 2 := (S 2 ) - {S} 2 « (S 2 ) (3.61) 



/oo roo roo roo 

df / df / dk / 
-oo J —oo J —oo J —oo 



dk' 



x*r(/ - f')5 T (k - k') (nl(f)n 2 (f)nl(k)h 2 (k')) Q(f)Q(k') (3.62) 



/oo roo roo rc 

df / df / dk / 
-oo J —oo J —OO J—' 



dk' 



x5 T (f - f')5 T (k - k') (^(/^(-^(^(-/^(AO) Q(f')Q(k') , (3.63) 

where we used the statistical independence and reality of ni(t) and n 2 {t) to obtain the last 
line. 

In Sec. |n A| , we defined the noise power spectrum Pf(|/|) in terms of the expectation 
value (rii(t)ni(f)) of the time domain random variables n,(t). (See Eq. ( |3.23|) .) An analogous 
expression holds in the frequency domain as well. Using definition (|3.49| ) for the Fourier 
transform hi(f) and definition ( 3.23Q for the noise power spectrum Pj(|/|), it follows that 



(nXf)W)) = \Kf ~ f) H\f\) ■ (3-64) 



Substituting this result into fl3.63p then yields 

1 roo roo 

a 2 ^-\ df df 6Uf - ma/i wi/'i)W) wo ( 3 - 65 ) 

4 J —oo J —oo 
T roo 

*jj $Px{\f\m\f\)\Q{f)W (3-66) 



where we replaced one of the finite-time delta functions 5x{f — f) by an ordinary Dirac 
delta function, and evaluated the other at / = /' to obtain the last line. 



16 For example, suppose we measure the outputs of the two detectors for a total observation time 
of one year. To estimate the cross-correlation (si(t)s2(t')), we simply form the products of si(t) 
and S2(t') for all t and t' having the same At = t — t' (e.g., 1 msec), and then average the results. 
We then repeat this procedure for At = 2 msec, 3 msec, etc. 



21 



To summarize: 



^ = IS T I df ^(1/1)7(1/1) W) , (3-67) 
T r°° 

a 2 ^- df P 1 (\f\)P 2 (\f\)\Q(f)\ 2 . (3.68) 

4 J — oo 

The problem now is to find the filter function Q(f) that maximizes the signal-to- noise 
ratio ( |3.51|) , with \i and a 2 as given above. This turns out to be remarkably simple if we 
first introduce an inner "product (A,B) for any pair of complex functions A(f) and B(f). 
The inner product of A(f) and B(f) is a complex number defined by 

/oo 
df • (3-69) 

-oo 

Since Pi(\f\) > 0, it follows that (A, A) > 0, and (A, A) = if and only if A(f) = 0. In 
addition, (A, B) = (B, A)* and (A, B + AC) = (A, B) + X(A, C) for any complex number 
A. Thus, (A, B) is a positive- definite inner product. It satisfies all of the properties of an 
ordinary dot product of vectors in three-dimensional Euclidean space. 

In terms of this inner product, the mean value /i and variance a 2 can be written as 

^-20^ r'i/i 3 A(i/i)p 2 (i/i)j [6J[}) 

a 2 ^^(Q,Q). (3.71) 

The problem is to choose Q(f) so that it maximizes the signal-to-noise ratio ( |3.51| ), or 
equivalently, the squared signal-to-noise ratio 

snr 2 = 4 » I ) t . ( 3 . 72 ) 

^ 2 (Q,Q) 

But this is trivial! For suppose we are given a fixed three-dimensional vector A, and are 
asked to find the three-dimensional vector Q that maximizes the ratio (Q-A) 2 /Q-Q. bmce 
this ratio is proportional to the squared cosine of the angle between the two vectors, it is 
maximized by choosing Q to point in the same direction as A. The problem of maximizing 
( P-72j ) is identical. The solution is 

n(f\-y 7(l/l)flgw(|/|) f o 7 oN 
WJ " l/l^d/IWI/l) ' [ } 

where A is a (real) overall normalization constant. 

One of the curious things about expression Q3.73j ) for the optimal filter Q(f) is that 
it depends upon the spectrum fi gw (/) of the stochastic gravity- wave background. This is 
a function that we do not know a priori. In practice this means that we cannot use a 




17 See, however, the discussion surrounding Eqs. (3.5S|) and (3.5 
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single optimal filter when performing the data analysis; we will need to use a set of such 
filters. For example, within the bandwidth of interest for the ground-based interferometers, 
it is reasonable to assume that the spectrum is given by a power-law fi gw (/) = Q a f a 
(where fl a = constant). We could then construct a set of optimal filters Q a (f) ( sa y, f° r 
a = —4, —7/2, ■ ■ • , 7/2, 4) with the overall normalization constants X a chosen so that 

H = tt a T . (3.74) 

With this choice of normalization, the optimal filter functions Q a {f) are completely specified 
by the exponent a, the overlap reduction function, and the noise power spectra of the two 
detectors. We would then analyze the outputs of the two detectors for each of these filters 
separately. Figure |] shows the optimal filter functions (displayed in the frequency domain) 
for both the initial and advanced LIGO detector pairs, for a stochastic background having 
a constant frequency spectrum O gw (/) = Q (i.e., a = 0)Q Figure ^| shows these same 
optimal filter functions displayed in the time domain 1.6., cLS cL function of the lag t — t! . 



Optimal Filter Functions 

(displayed in the frequency domain) 
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FIG. 4. Optimal filter functions Q(f) for the initial and advanced LIGO detector pairs, for 
a stochastic background having a constant frequency spectrum Q gw (/) = f2 . Both filters are 
normalized to have maximum magnitude equal to unity. 



Having found the optimal choice of filter function Q(f), it is now straightforward to 
calculate the signal-to-noise ratio for a given pair of detectors. Substituting ( |3.73D into 
( p.72j ) and taking the square root gives 
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The a in £l a and f a is just a number; it is not an index label. 
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gures f2l}]30 in Sec. [VI C| show the analogous optimal filter functions for different detector 
pairs. 
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Optimal Filter Functions 

(displayed in the time domain) 
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FIG. 5. Optimal filter functions Q(t — t') for the initial and advanced LIGO detector pairs, 
for a stochastic background having a constant frequency spectrum f2 gw (/) = Qq. Both filters are 
normalized to have maximum magnitude equal to unity. 
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T 
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f 6 Pl(\f\)P2(\f\) 



1/2 



(3.75) 



We will use this result in later sections to calculate signal-to-noise ratios and sensitivity 
levels for different detector pairs, assuming that the stochastic gravity-wave background has 
a constant frequency spectrum fi gw (/) = Q Q . Tables |I|-|V| in Sec. [VI D| contain the results of 
these calculations. 



IV. DETECTION, ESTIMATION, AND SENSITIVITY LEVELS 

Once the detectors have gone "on-line" and are generating data that needs to be an- 
alyzed, we will be confronted with the following questions: (i) "How do we decide, from 
the experimental data, if we've detected a stochastic gravity-wave signal?" (ii) "Assuming 
that a stochastic gravity- wave signal is present, how do we estimate its strength?" (iii) 
"Assuming that a stochastic gravity-wave signal is present, what is the minimum value of 
fin required to detect it 95% of the time?" In this section, we answer these questions, using 
a frequentist approach to the theory of probability and statistics-^ 



20 There are actually two approaches that one can take when analyzing data: the frequentist (or 
classical frequency probability) approach, which is adopted in this paper, and the Bayesian (or 
subjective probability) approach, which is adopted in Refs. [18,11]]. Although we will not describe 
the similarities and differences of these two approaches in any detail in this paper, it is important 



24 



A. Statistical considerations 



When performing a search for a stochastic background of gravitational radiation, it is 
convenient to break the data set (which might be hours, days, or weeks in length) into 
shorter stretches of N points, which we can then FFT and correlate with data from other 
detectors. Depending on the choice of N and on the sampling rate of the detector, these 
shorter segments of data will typically last on the order of seconds.^ From Eq. ( 3.50 ), we 
see that in a measurement over a single observation period T w 4 sec, the signal S is a sum 
(over / and /') of approximately 400 statistically independent random variables (products 
of the Fourier amplitudes of the signals). This is because Si(/) and s 2 (f) are correlated 
only when \ f — f'\ < 1/T « .25 Hz, and the bandwidth over which the integral in ( |3.50| ) gets 
its major contribution is ~ 100 Hz wide. Thus, by virtue of the central limit theorem, S is 
well-approximated by a Gaussian random variable, provided we are not too far away from 
the mode of the distribution. Equivalently, the values of S in a set of measurements over 
statistically independent time intervals (each of length T) are normally distributed. The 
mean value of this distribution is \x := (S), and the variance is a 2 := (S 2 ) — {S) 2 . 

Let s := (Si, S 2 , ■ ■ ■ , S n ) be a set of such measurements over statistically independent 
time intervals, each of length T.0 We can think of these measurements as n independent 
samples drawn from a normal distribution having mean \i and variance a 2 . The set s 
represents the outcome of a single experiment. From these samples, we can construct the 
sample mean 

/i:=if> (4.1) 

" i=i 



and the sample variance 



2 



a 



1 n 

— E(S-A) 2 . ( 4 - 2 ) 



n - i=i 



to emphasize that the frequentist and Bayesian approaches are inequivalent methods of analyzing 
data. Frequentists and Bayesians ask different questions about data and hypotheses, and conse- 
quently obtain different answers and draw different conclusions. In fact, there are certain questions 
that one can ask and answer in the Bayesian approach, that are ill-defined for a frequentist. In- 



terested readers should see Refs. [|l£4£] for a detailed discussion of the frequentist and Bayesian 



approaches applied to gravitational-wave data analysis with multiple detectors. 

21 For example, for N = 65536 = 2 16 and a sampling rate of 20 kHz, T = 3.2768 sec. For the 
same N and a sampling rate of 16.384 kHz = 2 14 Hz, T = 4.0 sec. Since the relevant bandwith for 
stochastic background detection is / < 300 Hz, the data stream could, in principle, be decimated 
to sampling rates which are substantially smaller — i.e., 1024 Hz. 

22 In order that the measurements be statistically independent, the time intervals should be non- 
overlapping, and T should be 3> the light travel time d/c between the two detectors. For the LIGO 
detector pair, this corresponds to T> 10~ 2 sec, which is satisfied for T w 4 sec. 
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Given the values of these estimators p, and a 2 , we would like to decide, in some reliable way, 
whether or not we have detected a stochastic gravity-wave signal. 

To make such a decision, we will apply a standard theorem from the theory of of proba- 
bility and statistics (see, e.g., p. 178 of Ref. pQ|). The theorem states that if fi is the sample 
mean of a set of n independent samples drawn from a normal distribution having mean fi 
and variance a 2 , then 

* = £f£ (4.3) 
o 7\/n 



is the value of a random variable having Student's t-distribution with parameter v = n — 
1. This is the classic Student's t-test. It is used to compare the means of two normal 
distributions that have the same variance, when the variance a 2 is unknown. 

Since tables of Student's t-distribution and its associated cumulative probability distri- 
bution function can be found in most handbooks on statistics (see, e.g., we could do 
all of the remaining calculations in this section in terms of the t-distribution. The drawback 
to this approach, however, is that the t-distribution depends on the parameter v — n — 1. 
This means that all of our results would depend on the number of observations n that con- 
stitute a single experiment. For stochastic background searches, this undesirable feature 
can be avoided by choosing n large enough so that the t-distribution and standard normal 
distribution are virtually indistinguishable. Since a typical total observation time will be on 
the order of months or years (~ 10 7 seconds), while T (the duration of a single observation 
period) is typically on the order of seconds, it is no problem to choose n ~ 10 3 or more.0 
For such large n, Eq. (|4.3|) can be rewritten as 

a/y/n 

where z is the value of a random variable having the standard normal distribution — i.e., z is a 
Gaussian random variable having zero mean and unit variance. Note that the approximation 
becomes a strict equality if a is replaced by its expected value a. This is because a linear 
combination of n Gaussian random variables (e.g., fi) is also a Gaussian random variable, 
independent of n. 

In the calculations that follow, we will want to assign probabilities to different events. 
From a frequentist point of view, this means that we should perform (or imagine performing) 
some fixed experiment many, many times. The probability of an event is then defined as 
the frequency of occurence of that event, in the limit of an infinite number of repeated, 
independent experiments. 



23 In fact, even for n = 30, the 2a values for the t-distribution and standard normal distribution 
differ by less than 5%. As n — > do, this difference goes to zero. 



26 



B. Signal detection 



In order to decide whether or not we've detected a stochastic gravity-wave signal, we 
need a rule that, given a set of measured data, will select for us one of two alternative 
hypotheses: 

H : A stochastic gravity- wave signal is absent. 

H\\ A stochastic gravity- wave signal is present, characterized by some fixed, but unknown, 
mean value \x > 0. 

Moreover, we would like this rule to be "optimal" with respect to some chosen set of criteria. 
This method of decision making, or "hypothesis testing" as it is more formally called, is a 
well-studied branch of frequentist statistics. As such, we will not go into any of the details 
here. After making the appropriate definitions, we will simply state the rule that we adopt 
for our stochastic background searches and explain in what sense it is optimal. Interested 
readers should see Ref. |2^| for a much more thorough discussion of the statistical theory of 
signal detection. 

To begin, let us note that the two hypotheses Hq and H\ are exhaustive and mutually 
exclusive — i.e., a stochastic gravity-wave signal is either absent or present. And, if present, it 
will be characterized by some fixed mean value /i, which is proportional to Qq for a stochastic 
background of gravitational radiation having a constant frequency spectrum £)gw(/) = ^o- 
The only requirement is that fi > 0. Hq is a simple hypothesis, since it does not depend on 
any unknown parameters. H\ is a complex (or composite) hypothesis, since it depends on a 
range of the unknown parameter \i. Explicitly, 



#i=U H, 



ndfj,,dfj, ; 



(4.5) 



n=0 



where is the hypothesis that a stochastic gravity- wave signal is present, characterized 

by a fixed mean value lying in the range (/i, \i + d/i] . 

As before, let s := (Si, S2, ■ • ■ , S n ) be a set of n statistically independent measurements 
of the cross-correlation signal S. Due to the noise intrinsic to the detectors and to errors 
inherent to the measurement process, the outcome of an experiment s is a random variable. 
It is described statistically by the probability density functions: 

p(s|0): Probability density function for the outcome of an experiment to be s, given that a 
stochastic gravity-wave signal is absent. 

p(s|/i): Probability density function for the outcome of an experiment to be s, given that a 
stochastic gravity-wave signal is present, characterized by the fixed mean value \l > 0. 



From (4.4), it follows that 

p(s|0) = (27r<T 2 )~ n/2 exp 

p{s\fi) = {2na 2 )- n/2 exp 



n q2 

k^ 2 \ 



E 

i=l 



(Si - 



To 2 



(4.6) 
(4.7) 
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The fact that we can use the sample variance a 2 (instead of the true variance a 2 ) on the 
right-hand-sides of (|4.6| ) and ( [4.7|) follows from the large n approximation that we used to 
obtain Eq. qPQ . 

A decision rule that, given the outcome of an experiment, selects for us either H or 
Hi is equivalent to a division of the space of all possible experimental outcomes into two 
disjoint regions Rq and R\\ If s 6 Ro, then Ho is chosen; if s £ Ri, then Hi is chosen. The 
success and failure of such a rule is characterized by two types of errors: A type I (or false 
alarm) error occurs when the decision rule chooses Hi when H is really true. A type II (or 
false dismissal) error occurs when the decision rule chooses H when Hi is really true. In 
terms of p(s\0), p(s\fi), R , and Ri, we have 

a = False alarm rate := / ds p(s\Q) (4-8) 

Jr.! 

/5(/i) = False dismissal rate := / ds p(s\fi) (4.9) 

JRo 

Note that, for the complex hypothesis Hi, the false dismissal rate is actually a function of 
the mean value fi > 0. Note also that 1 — a is the fraction of experimental outcomes that the 
decision rule correctly identifies the absence (not presence!!) of a stochastic gravity- wave 
signal. If we want to talk about detection, we should evaluate 

7(/i) = Detection rate := 1 — /3(/i) , (4-10) 

which is the fraction of experimental outcomes that the decision rule correctly identifies the 
presence of a stochastic gravity-wave signal, characterized by the fixed mean value fi > 0.Q 
In order for the decision rule to choose the regions Rq and Ri in an "optimal" way, 
we must first select some set of criteria with respect to which "optimal" can be defined. 
For stochastic background searches, where one does not know a priori the "costs" that one 
should associate with false alarm and false dismissal errors, it is reasonable to choose a 
decision rule that minimizes the false dismissal rate /3(yu) for a fixed value of the false alarm 
rate a. Equivalently, one chooses a decision rule that maximizes the probability of detecting 
a stochastic gravity-wave signal, while keeping the false alarm rate fixed. This decision 
criterion is known in the literature as the Neyman-Pearson criterion. In general, for a 
complex hypothesis, the Neyman-Pearson criterion yields regions Ro and Ri that depend on 
the unknown parameter /i. But for this case, where Hi involves all parameter values fi > 0, 
Ro and Ri are actually independent of Without going into details here, let us simply 
state the result: Namely, in the context of stochastic background searches as described 
above, the Neyman-Pearson criterion is satisfied if, given the outcome of an experiment s, 
we form the estimators fi and a 2 according to Eqs. ( |4.1| ) and ( |4.2| ), and then 



choose Ho if fi < z a fr/^/n, and choose Hi if fi > z a a / ' -Jn . (4-11) 



24 The detection rate 'y(n) should not be confused with the overlap reduction function 7(/), which 
was defined in Sec. |III B| . 

25 See pp.152-153 of Ref. p| for more details. 
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Here z a is that value of the random variable z for which the area under the standard normal 
distribution to its right is equal to a. (See Fig. ^.) In terms of the complementary error 

Standard Normal Distribution 

0.5 
0.4 
0.3 
0.2 
0.1 



-4.0 0.0 z 4.0 

z 

FIG. 6. The standard normal probability distribution for a Gaussian random variable having 
zero mean and unit variance. z a is that value of z for which the area under the standard normal 
distribution to its right is equal to a. Typical values for 1 — a are 0.90, 0.95, and 0.99. The 
corresponding values for z a are zn.io = 1-28, zn.05 = 1.65, and zn.oi = 2.33. These are the threshold 
values appropriate for the one-sided test described in the text. 

function 

erfc(z) := — / dx e~ x (4.12) 

y/Tl Jz 

(see Fig. 0), 

z a = V2 erfc~ 1 (2a) , (4.13) 

where erfc -1 denotes the inverse of erfc. Since fx and a are functions of the outcome of the 
experiment s, the inequalities fi < z a a/y/n and fi> z a a/y/n define Rq and R\. Note also 
that the decision rule can be restated as 

choose Ho if y/n SNR < z a , and choose Hi if y/n SNR > z a , (4-14) 

where 

yfn~ SNR := v 7 ^ - (4.15) 
a 

is the measured signal-to-noise ratio after n observation periods. The fact that the signal- 
to-noise ratio enters the above inequalities is one of the main reasons why we paid so much 
attention to evaluating it in Sec. [TTT|. 




29 



Complementary Error Function 



FIG. 7. The complementary error function. 



Given the above decision rule for signal detection, we can now calculate the false alarm 
and false dismissal rates defined by Eqs. ( j4.8|) and (|4.9| ). First, for the false alarm rate, one 
finds 

a = Prob(/i > z a &/y/n\n = 0) (4.16) 
= Piob(z>z a ) (4.17) 




(4.18) 

where we used Eq. ( [4.4|) with \i — to obtain the second equality. Thus, our two uses of the 
symbol a (for the false alarm rate and for the area under the standard normal distribution 
to the right of z a ) are consistent. 

Second, for the false dismissal rate, 

P(jjt) = Prob(/i < z a ojyfn\\i > fixed) (4.19) 
= Prob(2; < z a - y/n fi/a) (4.20) 




V2 



(4.21) 



where we again used Eq. (|4.4|) (but this time with /i > 0) to obtain the second equality. A 
graph of the detection rate y(p) '■— 1 — is shown in Fig. |8|. From this graph, we see 
that the detection rate 7(/x) approaches a as ji — » 0; it approaches 1 as ii — > oo. Also, the 
detection rate equals 0.50 for \i = z a cr/y/n (or, equivalently, for y/n fi/fr = z a ). Thus, if we 
replace a with its expected value a, we see that the detection rate is only 50% for a signal 
having a theoretical signal-to-noise ratio after n observation periods equal to z a . 

It is also interesting to note that the Neyman- Pearson detection criterion, when applied to 
stochastic background searches, is equivalent to the maximum-likelihood detection criterion 
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Detection Rate 



1.0 




FIG. 8. The detection rate 7(//) := 1 — /?(//) plotted as a function of the mean value /i. 



In other words, we will obtain the same decision rule given above if we first construct 
the likelihood ratio 

p(s\fi) 



a(s|a*) 



p(s|0) 



(4.22) 



(4.23) 



maximize A(s|/i) with respect to variations of the parameter /i > 0: 

A max (s) := max A(s|/x) ; 

£t>0 

and then divide the space of all possible experimental outcomes s into two regions Rq and 
Ri according the the rule 

choose H if A max (s) < An, and choose Hi if A max (s) > A , (4.24) 

where An is chosen so that the false alarm rate equals a. Explicitly, from Eqs. (|4.6|) and 
(|3), it follows that 



A(s|/x) = exp 



exp 



na 

^2-1/' -/' 



(4.25) 



A(s|/i) is maximized when /i equals the sample mean /t:Q 

1 nfi 2 



A max (s) = A(s|/t) = exp 



2 a 2 



(4.26) 



The decision surface A is obtained by setting /t equal to the threshold value z a a/y/n: 



An = A 



maxima 



a/y/n) = exp 



(4.27) 



26 The sample mean /t is said to be the maximum-likelihood estimator for this problem. 
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C. Parameter estimation 



Assuming that a stochastic gravity- wave signal is present, characterized by some fixed, 
but unknown, mean value \i > 0, parameter estimation attempts to answer the question: 
"What is the value of /i?" It does this by first constructing the sample mean fi and sample 
variance a 2 of a set s := [Si, S 2 , ■ ■ ■ , S n ) of n statistically independent measurements of the 
cross-correlation signal S, as described in Sec. IV A . From Eq. (O), we then know that 



a 



n 



(4.28) 



is the value of a Gaussian random variable having zero mean and unit variance. Thus, if 
z a /2 is that value of the random variable z for which the area under the standard normal 
distribution to its right is equal to a/2 (see Fig. ||), then (1 — a) ■ 100% of the time 



Standard Normal Distribution 




FIG. 9. The standard normal probability distribution for a Gaussian random variable having 
zero mean and unit variance. z a / 2 is that value of z for which the area under the standard normal 
distribution to its right is equal to a/2. The area under the standard normal distribution between 
~ z a/2 an d z a /2 is thus 1 — a. Typical values for 1 — a are 0.90, 0.95, and 0.99. The corresponding 
values for z a / 2 are zn.05 = 1.65, zo.025 = 1-96, and £0.005 = 2.58. 



Z a/2 < Z < Z, 



a/2 , 



(4.29) 



or, equivalently, 



fi - z a/2 a I \/n < ji< fi + z a/2 a I y/n . 



(4.30) 



Said another way, in an ensemble of observations of the same stochastic background, a 
fraction 1 — a of the intervals 



I a :-- 



jj, - z a/2 (r/y/n, jj, + z a /2 a/y/n 



(4.31) 
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constructed from the measured data, will contain the value of the true mean fi. Equivalently, 
a is the fraction of intervals I a that fail to contain the value of the true mean \i. Of course, 
given the outcome of a single experiment s, the interval I a either contains or does not contain 
the value of /i. And the value of /z, if it is contained in I a , need not be any closer to the center 
of the interval than to either of its edges. Thus, the confidence that one associates with the 
above estimation procedure is not equivalent to our degree of belief that the true mean [i lies 
within a given interval. (This is a Bayesian interpretation of probability.) Rather, it is the 
fraction of experimental outcomes that our estimation procedure will produce an interval 
that contains the true mean /i, in the limit of an infinite number of repeated, independent 
experiments. 



D. Sensitivity levels 

Let us assume once again that a stochastic gravity- wave signal is present, characterized 
by some fixed, but unknown, mean value /i > 0. Then it is reasonable to ask: "What is the 
minimum value of fi required so that our decision rule correctly identifies the presence of a 
signal at least 7 • 100% of the time?" 

The answer to this question can be obtained by applying the results of Sec. |1V Bl Namely, 
we require that the detection rate 7(/x) := 1 — /3(/i) be greater than or equal to the desired 
rate 7, and then solve the resulting inequality for fi. Explicitly, from Eq. (|4.21 ), it follows 
that 

l-*) = ie, f c(^|^)> 7 , (4.32) 

or, equivalently, 

z a - y/E fijo < V2 erfc _1 (27) , (4.33) 

where we replaced the estimator a by its expected value a to obtain the LHS of the above 
equation. Thus, 

fi> (z a -V2 erfc- 1 (2 7 )) (4.34) 
= ^=V2 (erfc _1 (2a) - erfc" 1 ^)) . (4.35) 



Equivalently, 



v 7 ^ SNR > V2 (erfc" 1 (2a) - erfc" 1 (2 7 )) , (4.36) 



where y/n SNR := y/n fi/a is the theoretical signal-to-noise ratio after n observation periods. 
Note that the right-hand-sides of Eqs. ( |4.35| ) and ( |4.36| ) depend on both a and 7.Q This 



27 The dependence on a is via the threshold value z a , which defines the decision rule (4.14) 
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means that to calculate the minimum value of fi (or the minimum signal-to-noise ratio), 
we must specify the desired detection rate 7 in addition to the false alarm rate a. In the 
past (see, e.g., Refs. people have only specified the false alarm rate a. It seems that 

they have mistakenly assumed that the probability of correctly identifying the presence of a 
signal was 1 — a. But, as mentioned in Sec. [IV B| , 1 — a is the probability that the decision 
rule correctly identifies the absence of a signal. To talk about detection (and the associated 
minimum value of /1), one needs to include the additional parameter 7. 

For stochastic background searches, we can rewrite Eq. ( |4.35|) in terms of Qq, 7(/), and 
the noise power spectra P%(f) and of the two detectors: 



n > 



1 IOtt' 



df 



7 2 (|/|) 



PPl(\f\)W\) 



-1/2 



v / 2(erfc" 1 (2a) -erfc _1 (27)) , (4.37) 



where T tot := nT is the duration of a single experiment s. This result follows from Eq. (|3.75|) 
for a stochastic background of gravitational radiation having a constant frequency spectrum 
^gw(/) = ^o- Let us now evaluate the minimum value of Qq for 4 months of observation 
(i.e., T tot = 10 7 sec), for a false alarm rate a = 0.05, and for a detection rate 7 = 0.95. 
Denoting the solution by fl °, we fincQ 



(i) ^ 5% ' 5% hf 00 

(ii) fi"*** hl 00 



5.74 x 10 6 for the initial LIGO detector pair; 
5.68 x 10~ n for the advanced LIGO detector pair. 



Note, however, that these values disagree with those quoted in the literature (see, e.g., 
Refs J7|||). In Refs. the minimum value of Qq is determined by the equation 



lOvr 2 



df 



7 2 (l/l) 



/ 6 A(|/|)P 2 (|/|) 



'1/2 



V2 erfc _1 ( 



a) 



(4.38) 



There are two mistakes: First, the argument of the inverse complementary error function 
is a instead of 2a. But this is appropriate for a two-sided test, which would be correct 
if the mean value fi could be either positive or negative. Second, and more importantly, 
there is no term proportional to erfc - (27). As mentioned above, it seems that the authors 
mistakenly assumed that the probability of correctly identifying the presence of a signal was 
1 — a. Even the a = 0.10, 7 = 0.90 value of Eq. (|4.37|) does not agree with the values of 
fig 0% quoted in Refs. Thus, by Qq 0% the authors do not mean the minimum value of 

Qq for a false alarm rate equal to 10% and a detection rate equal to 90%. 

Alternatively, the absence of the term erfc -1 (27) from Eq. ([4.38|) is equivalent to calcu- 
lating the minimum value of Sl f° r a detection rate 7 = 0.50. This is because erfc _1 (l) = 0. 
Thus, the values of S7q 0% quoted in Refs. [7|,§ are for a false alarm rate equal to 5% and a 



28 There is nothing special about the choice a = 0.05 and 7 = 0.95. For example, a and 7 need not 
sum to one. We could equally well have chosen a = 0.10 and 7 = 0.95, and calculated a different 

■ ■ 1 O 95%,10% 

minimum value "0 
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detection rate equal to 50% — not 90% as they claim. This is why the minimum values of 
flo quoted in those papers are smaller than those found here. 

Table | in Sec. |V1D| contains theoretical signal-to-noise ratios after 4 months of obser- 
vation, for a stochastic background having a constant frequency spectrum fi gw (/) = Q = 
6 x 10 6 hToo- 

Table |lj in Sec. [VI D| contains minimum values of Qq ^ioo f° r 4 months of observation, 
for a false alarm rate equal to 5%, and for a detection rate equal to 95%. 



E. Summary 

We started this section by asking a series of questions. To conclude, we summarize the 
answers obtained above, and address a couple of other related issues. 

(i) "How do we decide, from the experimental data, if we've detected a stochastic gravity- 
wave signal?" 

Answer: We compare the measured signal-to-noise ratio after n observation periods to 
the threshold value z a . If sfn SNR < z a , we conclude that a stochastic gravity- wave signal 
is absent. If ■s/n SNR > z a , we conclude that a stochastic gravity-wave signal is present, 
characterized by some fixed, but unknown, mean value /i > 0. 

Note, however, that we can never conclude, with 100% confidence, that a stochastic 
gravity- wave signal is absent or present. Our decision rule leads us to infer one of these two 
possibilities, but the rule is not perfect. The false alarm rate a and false dismissal rates 
P(fi) (defined by Eqs. ( f4.8|) and (|4.9|) ) are the error rates associated with the rule. Thus, 
claims about the absence or presence of a stochastic gravity-wave signal should always be 
made with these error rates in mind. 

Moreover, one has to be very careful about trying to define "termination" criteria. For 
example, it would be misleading to try to terminate an experiment by correlating the out- 
puts of two gravity-wave detectors until the measured signal-to-noise ratio for the total 
observation period exceeds some threshold value z a . One can show that the false alarm rate 
associated with such a rule is 100%! In other words, the conclusion drawn from such an 
experiment would always be that a stochastic gravity- wave signal is present, even in the 
absence of a signal. Noise intrinsic to the detectors and errors inherent to the measure- 
ment process are sufficient to guarantee that the measured signal-to-noise ratio for the total 
observation period will eventually exceed any z a . 

If, however, the value of the threshold level increases with observation time in an appro- 
priate manner, then one can define termination criteria that have false alarm rates less than 
100%. A famous theorem of probability and statistics, called The Law of the Iterated Loga- 



rithm [23], states that if Si, S2, • • • are statistically independent and identically distributed 



random variables with zero mean and finite variance a 2 , thenP] 



29 The lim sup n ^ oc of a sequence xi, X2, • • • is defined as follows: Let a m equal the least upper bound 
of the subsequence x m ,x m+ i, ■ ■ ■ . (Note that ai > 0,2 > • • • •) Then limsup,^^ x n := lmim^oo a 



ni ■ 
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Pl ,,l, [ limsup Si + S 2 + --- + Sn = ! | = j ( | 

This means that if A < 1, there is unit probability that 

S i + S * + --- + S - > A (440) 
ay/2n log(log(n)) 

for infinitely many n. If A > 1, there is unit probability that (|4.40 ) holds for only finitely 



many n. In terms of the measured signal-to-noise ratio after n observation periods, Eq. ( |4.40 ) 
can be rewritten as 



SNR > A ^/21og(log(n)) , (4.41) 

where we have replaced a in the definition of SNR by its expected value a. Thus, if we want 
to define a termination criterion with a false alarm rate less than 100%, we should compare 
/ n SNR with threshold levels t n satisfying 

lim tn = oo . (4.42) 
/log(log(n)) 



This is consistent with the claim made in the previous paragraph that constant threshold 
levels t n = z a will always be exceeded for some n. Unfortunately, we have not been able 
to write down a simple analytic expression for the false alarm rate, for arbitrary threshold 
levels t n satisfying (|4.42|) . The probability that y/n SNR > t n for some n involves integrals 
of products of the Gaussian probability density function with the complementary error 
function. 

(ii) "Assuming that a stochastic gravity-wave signal is present, how do we estimate its 
strength?" 

Answer: Assuming that a stochastic gravity- wave signal is present, characterized by some 
fixed, but unknown, mean value /x > 0, we estimate /i by constructing the interval 



a • 



fi - z a/2 oj y/n, ft + z a /2 oj yfri . (4.43) 



In an ensemble of observations of the same stochastic background, 1 — a is the fraction of 
intervals constructed in this way that contain the value of the true mean /i. 

We should emphasize that parameter estimation assumes the presence of a signal. First, 
as mentioned above, we can never be 100% certain that a stochastic gravity-wave signal is 
present. Second, it does not make any sense to try to estimate the parameters of something 
that we assume does not exist. 

(iii) "Assuming that a stochastic gravity-wave signal is present, what is the minimum value 
of Qq required to detect it 7 ■ 100% of time?" 
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Answer: Assuming that a stochastic gravity- wave signal is present, characterized by some 
fixed, but unknown, mean value /i > 0, the minimum value of Qq, for an observation time 
T tot , false alarm rate a, and detection rate 7, is given by 



n 



1 



IOtt 2 



\Mtot 3-f/o 



df 



7 2 (|/|) 



/ 6 A(|/|)P 2 (|/|) 



1 2 

For fixed a and 7, the factor of T tot implies that the minimum value of £Iq decreases with 
increasing observation time. Thus, the sensitivity improves as the total observation time 
increases. This means that as T tot increases we can put a tighter "upper limit" on the values 
of Qq for stochastic gravity- wave signals that we will falsely dismiss more than (1 — 7) • 100% 
of the time. 

1/0 

Alternatively, for a fixed Qq and a fixed false alarm rate a, the factor of T tot implies that 
the false dismissmal rate 1 — 7, for a stochastic gravity-wave signal having a strength equal 
to Qq, decreases with increasing observation time. For example, suppose that a particular 
theory — like cosmic strings — predicts a value of Qq = 10~ 7 . Moreover, suppose that, in 
successive years of observation, we fail to detect the presence of a signal at this level of 
sensitivity. Then we can still say that the probability of our falsely dismissing a stochastic 
gravity-wave signal having a strength equal to flo = 10~ 7 has decreased over the course of 
the observation. 

The above upper limit on stochastic background signal strengths (defined by Eq. ( |4.44Q ) 
is different from that obtained by setting the measured signal-to-noise ratio for the total 
observation period equal to the threshold value z a : 



-1/2 



v / 2(erfc" 1 (2a) -erfc" 1 (27)) . (4.44) 



1 



tot 



IOtt 2 

3m 



df 



7 2 (l/|) 



f 6 Pi(\f\)P2(\f\) 



-1/2 



V2 erfc" 1 (2a) 



(4.45) 



This alternative upper limit is the maximum value of the test statistic leading to the con- 
clusion that a stochastic gravity- wave signal is absent. The decision rule — by its very 
construction — will never allow one to conclude that a stochastic gravity-wave signal is present 
with an Qq less than this upper limit. 

The above two definitions of upper limit agree when the detection rate 7 = 0.50. This is 
because erfc~ 1 (l) = 0. This means that, if a stochastic gravity-wave signal is present, with an 
Qo less than or equal to this maximum value, then there is more than a 50% chance of falsely 
dismissing it. Note also that if we change the argument of the inverse complementary error 
function from a to 2a in Eq. (|4.38|) , then Eqs. ( |4.38| ) and ( |4.45|) agree. Thus, the minimum 
values of flo quoted in Refs. |7]|| can be interpreted as either the 7 = 0.50, a = 0.05 upper 
limit defined by Eq. ( |4.44j ), or the a = 0.05 upper limit defined by Eq. ( |4.45|) . 



V. COMPLICATIONS 



In Sec. m, we discussed optimal filtering under the assumptions that the noise intrinsic 
to the detectors were: (i) stationary, (ii) Gaussian, (iii) statistically independent of one 
another and of the stochastic gravity-wave background, and (iv) much larger in magnitude 
than the stochastic gravity-wave background. In this section, we describe the modifications 
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that are necessary when some of these assumptions are removedf°| We will also describe 
how one can correlate the outputs of 4 or more detectors, and how one can combine data 
from multiple detector pairs to increase the sensitivity of a stochastic background search. 



A. Signal-to-noise ratios for arbitrarily large stochastic backgrounds 

In Sec. [Ill (J| , we calculated the signal-to-noise ratio for a stochastic background of grav- 
itational radiation, assuming that the noise intrinsic to the detectors were much larger in 
magnitude than the gravitational strains. Although this assumption is most likely valid, 
there are at least two reasons why we want to consider stochastic background signals whose 
magnitudes are comparable to (or larger than) the noise intrinsic to the detectors: (i) Com- 
puter simulations for stochastic background searches allow one to "dial-in" arbitrarily large 
signal strengths. Thus, in order to compare theoretical predictions with the results of com- 
puter simulations, we need to be able to analyze the large signal case, (ii) Future generations 
of gravity-wave detectors might have intrinsic detector noise levels comparable to the level 
of a real stochastic background signal. For this case, optimal filtering of large signal data 
will also be necessary. 

To begin, let us recall the main results of Sec. Q Under the assumption that the noise 
intrinsic to the detectors were much larger than the gravitational strains, the calculation of 
the variance a 2 := (S 2 ) — (S) 2 simplified considerably. For the large noise case, a 2 was 
dominated by the pure detector noise contribution: 

^7/ #^i(l/l)Wl)lW)l a - (5.1) 

4 J—oo 



This is Eq. (ETBSQ from Sec. PTTC 



If, however, the magnitude of the stochastic background signal is comparable to the noise 
intrinsic to the detectors, the calculation of the signal-to-noise ratio SNR := fi/a is more 
involved. Although the mean value \i is independent of the relative size of the stochastic 
background signal and the detector noisejf] 

f=^ T IJf fi S w(|/|) 7(1/1) W) , (5-2) 



the variance a is not. Explicitly, 



4 

where 



J 1 roo 

o 2 = -r\ df \Q(f)\ 2 R(f) , (5.3) 



30 We will always assume that the noise intrinsic to a detector is Gaussian and statistically inde- 
pendent of the gravitational strains. 
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See the discussion surrounding Eqs. ( 3.52| )-( |3757 ) in Sec. Ill C . 
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R(f) := Pi(|/|)P 2 (|/|) + TTT^ ~ JE S^ (Ad/I) + ^(1/1) 




(5.4) 

Note the additional terms that contribute to the variance. Roughly speaking, they can be 
thought of as two "signal+noise" cross-terms and a "pure signal" variance term.0 When 
^gw(|/|) is negligible compared to the detector noise power spectra Pi(\f\) (i = 1,2), 
Eqs. ( |5.3j ) and ( |5.4j ) for a 2 reduce to the "pure noise" variance term ( |5.1| ) as they should. 

At this stage of the analysis, the filter function Q(f) in (|5.2|) and ( |5.3|) is arbitrary. For 
the case of large detector noise, we were able to make an optimal choice for Q(f), which 
maximized the signal-to-noise ratio. This was facilitated by introducing an inner product 

/oo 
df A*(/W)Pi(|/i)^(m) , (5.5) 
-oo 

and writing \i and a 2 in terms of this inner product. The squared signal-to-noise ratio then 
took the form 

SNR 2 = ( Ml ) T ^ , (5.6) 

which was maximized by choosing 

Although Q(f) depended on f2 gw (/), we could construct a set of optimal filters Q a (f) for 
stochastic backgrounds having power-law spectra Q^f) = VL a f a (where Q a = constant). 
The proportionality constants Q a could always be absorbed into the normalization constants 
X a . The resulting set of optimal filters were then completely specified by the exponent a, 
the overlap reduction function, and the noise power spectra of the two detectors. 

For the case where the stochastic background is comparable to the noise intrinsic to the 
detectors, we can try to do something similar. We can define a new inner product 

/oo 
df A*(f)B(f)R(f) , (5.8) 
-oo 

where R(f) is given by fl5.4Q . Although this inner product is more complicated in form than 
the original inner product ( |5.5| ), it is still positive-definite, since R(f) is real and positive. 
In terms of (|5.8|) , 

^=^T(Q,A) and a 2 = |(g,Q), (5.9) 



32 These are the terms proportional to f2 gw (|/|) and f2;L.(|/|), respectively. 
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where 



A(f) ■■= i/r 3 ^(\f\H\f\)R-\f) 

The squared signal-to-noise ratio is thus 

{Q,Af 



SNR 



m^ 2 



IOtt 2 



T 



(Q,Q) 



(5.10) 



(5.11) 



which has the same form as ( |5.6|) , although with a much more complicated expression for 
A(f). But the same argument for maximizing the squared signal-to-noise ratio still goes 
through. The optimal choice of filter function is 



Q(f) = XA(f) 



(5.12) 



where A is a (real) overall normalization constant. 

But this is where the similarity with the large detector noise case ends, and where the 
complications start to arise. The main problem is that the optimal filter function Q(f) has 
a complicated functional dependence on the stochastic gravity- wave spectrum fi gw (/). For 
the large detector noise case, this problem did not exist. As mentioned earlier, we could 
always consider stochastic backgrounds having power-law spectra Qgw(f) = ^af a , and then 
construct a set of optimal filters Q a (f) labeled by the different values of a. But for the 
optimal filter function Q(f) given by ( |5.12| ), ( |5.10| ), and (|5.4|) , this idea of constructing a 
set of filters labeled by only the power-law exponents a fails. The proportionality constants 
Q a cannot be absorbed into the normalization constants A a . For example, if we consider 
a stochastic background having a constant frequency spectrum Q gw (f) = Q (i.e., a = 0), 
then 



Q(f) = Afiol/r 7(1/1) 



IOtt* \f\* 



Pi{\f\) + P*(\f\) 



+ 



3HI 
IOtt 2 



f (i + r(l/l) 



(5.13) 



Although the factor of Qq in the numerator can be absorbed into the normalization constant 
A, the factors of Qo an d m the denominator cannot. In other words, if we want to construct 
a set of optimal filters for arbitrarily large stochastic backgrounds having power-law spectra 
^gw(/) = ^af a , we need to specify the proportionality constants fl a in addition to the 
exponents a. The "space of optimal filters" thus becomes a much larger set, parametrized 
by (a,Q a ). Although in principle this poses no problem, in practice it requires a more 
sophisticated search algorithm; the detector outputs will have to be analyzed for each of the 
filters Q( a ,Q a )(f) separately. 



B. Nonstationary detector noise 

It is not at all uncommon for the power spectra of the noise intrinsic to the detectors to 
change over the course of time. There will be periods of time when the detectors are rela- 
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tively "quiet," and other periods of time when the detectors are relatively "noisy. "0 These 
changes in the power spectra will, in turn, lead to measurements whose statistical properties 
also change with time. For example, during the quiet periods of detector operation, the 
measurements ^'Si, ^S^, • • • , ^S ni will have an associated variance ^V 2 that will be smaller 
than the variance ( 2 V 2 associated with the measurements ^Si, ( 2 v?2, • • • , ^S n2 taken during 
the noisy periods. Moreover, if the optimal filter function is normalized (by an appropriate 
choice of A in Eq. ( |3.7$ )) so that the theoretical mean /i is equal to some fixed value (e.g., 
\i = fl a T for a stochastic background having a power-law spectrum O gw (/) = fl a f a ), then 
the quiet periods will have a correspondingly larger signal-to-noise ratio. Thus, a natural 
question that arises in this context is: "Can one combine the different sets of measurements, 
corresponding to the quiet and noisy periods of detector operation, so as to maximize the 
overall signal-to-noise ratio?" The answer to this question is "yes," and the proof is sketched 
below. 

In order to handle the most general case, let us consider m different sets of measurements 

«S 2 , ••• , ®S ni (5.14) 

corresponding to m different levels of detector noise, or m different periods of detector 
operation. (Here i = 1,2,- -^m.) Each of these measurements is taken over an identical 
time interval of length T. These measurements can be thought of as realizations of m random 
variables ^S, each having the same theoretical mean 

% := («£) =: Li , (5.15) 

but different theoretical variances 

(V := (^S 2 ) - (%) 2 . (5.16) 

The equality of the mean values follows because we assume identical normalization conven- 
tions on the optimal filters — e.g., li = Q a T for a power-law spectrum Q^^f) = Q a f a . We 
also assume that all of the measurements ^Si, ^HS^, • ■ ■ , (jn ^S nm are statistically independent 
of one another.^] 

For each set of measurements, we can construct the sample mean (or estimator) 



«a:=^£«V (5.17) 
I j=l 

Viewed as a random variable in its own right, W/i has mean value 



33 Although the frequency of these more "noisy" periods should decrease as the detectors are 
gradually improved over the course of months or years, variations in the detector noise power 
spectra will still inevitably occur. 

34 As mentioned in Sec. EVA, this means that the measurements are taken over distinct, non- 



overlapping periods of operation, with T 3> the light travel time d/c between the two detectors. 
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IH : = ((%) = fj, , (5.18) 

and variance 



2 



( ( % 2 ) " ( W A) 2 = — • (5-19) 



n, 



What we want to do now is combine all the measurements ^'Si, ^IS^, • ■ ■ , ^S nm (or, 
equivalently, combine the sample means ^fi, ^fi, ■ ■ • , ( m )/t) so as to maximize the overall 
signal-to-noise ratio. We thus use a weighted average to define the estimator 

:= V m X > ( 5 - 20 ) 

and then choose A, > to maximize the signal-to-noise ratio of fx. From ( 5.20| ), it follows 
that p, has mean value 

/^ := (A>=^ ( 5 - 21 ) 
(which is independent of the choice of Aj), and variance 

s-^m \2 ,_2 

4 := (A 2 ) - (A) 2 = ^ ? a • (5-22) 
The squared signal-to-noise ratio is thus 

SNRJ • (5.23) 

To find the Aj which maximize SNR| (and hence which tell us how to optimally combine 
data from periods of quiet and noisy detector operation) is quite easy. This is because 
Eq. (|5.23|) can be written as a ratio of inner products, just as we were able to write the 
squared signal-to-noise ratio ( p.72|) in Sec. [Ill Q as a ratio of inner products. Explicitly, 



SNR^ = ^ v ' ' (5.24) 
where the inner product (oc,(3) is defined by 

m 

(a, /?):=£ (5.25) 



i=i 



for any pair of (complex) sequences and /3j (i = 1,2, ■ • • ,m). The inner product (a, (3) 
is positive-definite, since of is real and positive, and it satisfies all of the properties of the 
ordinary dot product of vectors in three-dimensional Euclidean space. Thus, choosing 

Aj oc a," 2 (5.26) 
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— * — * 

maximizes ( |5.24p in the same way that choosing A proportional to B maximizes the ratio 
(A ■ B) 2 /(A ■ A). By averaging each set of measurements with the inverse of its associated 
theoretical variance, we give more weight to signal values that are measured when the 
detectors are quiet, than to signal values that are measured when the detectors are noisy. 
This weighting maximizes the overall signal-to-noise ratio. 

For the optimal choice of weights Aj = a~ 2 , the inverse variance of the optimal estimator 



^optimal AL= CT r 2 (5.27) 



has a very simple form: 



^optimal = E =£"i ( 'V 2 • (5-28) 

i=\ i=l 

This result says that the variances for the optimal estimator add like electrical resistors in 
parallel. The squared signal-to-noise ratio of the optimal estimator also has a very simple 
form: 

m m 

SNR 2 ptimal = £ SNR 2 = £ n t «SNR 2 , (5.29) 



i=l 



where 



and 



SNRj := ^ = ^ (5.30) 



WSNR:=^ = ^- . (5.31) 

Thus, the squared signal-to-noise ratio of the optimal estimator is simply a sum the squared 
signal-to-noise ratios for each after observation periods, each of length T. 

It is instructive to compare the results for the optimal estimator /2 op timai with those for 
the "naive" estimator 

Anaive := ^ "^3 ' (5.32) 

n tot j = i j = i 

which simply averages all of the signal estimates, paying no attention to the different vari- 
ances W<r 2 associated with the ^S. The naive estimator /t na i ve corresponds to Eq. (|5.20|) with 
Aj = rii and n tot := YhLi^- One can show that 

i m i m 

^naive = 4" £ O* = -5- £ <V (5.33) 

"tot j=l "tot j=l 



and 



1 m 1 m 

SNR- 2 ve = — £ "? SNRr 2 = — £^ W SNR- 2 . (5.34) 



n tot j=l n tot j=l 
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Note that in addition to the factors of l/n\ ot , the signs of the exponents for the variances 
and signal-to-noise ratios in Eqs. Q5.33 ) and ( 5.34Q are opposite to those in Eqs. ( 5.28|) and 

To give a numerical example, suppose we have just two different periods of detector 
operation, with the second period twice as noisy as the first— i.e., ( % = 2 «<7. Since 
of = ( l h 2 /ni, it follows that if there are four times as many measurements during the 
noisy period (i.e., n 2 = 4ni), then a\ = o~\. Intuition suggests that the optimal way of 
combining the measurements for this case is to weight the estimators ^fi and ^fi equally 
(i.e., Ai = A 2 = o~i 2 ). This agrees with the above mathematical analysis, and yields a 
squared signal-to-noise ratio for the optimal estimator equal to 

SNR 2 ptimal = 2n x «SNR 2 = 2m ^ . (5.35) 
For the naive estimator, 

SNR^ve = (1) SNR 2 = ^ . (5.36) 

Thus, for this particular example, the signal-to-noise ratio for the optimal estimator is 
1/34/25 = 1.17 times larger than that for the naive estimator. 



C. Multiple detector pairs 

Combining measurements from multiple detector pairs in order to increase the sensitivity 
of a stochastic background search is identical to combining different sets of measurements 
corresponding to quiet and noisy periods of detector operation in order to maximize the 
overall signal-to-noise ratio. The only difference between the two is one of interpretation 
and notation. In Sec. IV BL 



®Si, (i) S 2 , •■■ , % nj (5.37) 

denoted rij different measurements, each of length T, of the optimally-filtered cross- 
correlation signal ®S when the level of detector noise was characterized by the variance 

Waf : = («S 2 ) - («S> 2 . (5.38) 

In this section, for multiple detector pairs, 

«% 2 , ■•■ , (% ni . (5.39) 

denote riij different measurements, again each of length T, of the optimally-filtered cross- 
correlation signal W'S between the ith and jth detectors.^ As usual, 



In Sec. VB, i = 1, 2, • • • , m labeled m different levels of detector noise, or m different periods of 
detector operation. In this section, i,j = 1, 2, • • • , I label I different detectors, and m := 1(1 — l)/2 
is the number of different detector pairs. 
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denotes the variance of the cross-correlation signal W'S. 

Provided that each cross-correlation signal measurement occurs over the same time in- 
terval T (which should be 3> the light travel time between any pair of detectors), and that 
the optimal filter functions for each detector pair have identical normalizations (e.g., A in 
Eq. (|3.73|) is chosen so that = (^'S) = Vt a T for a stochastic gravity- wave background 
having a power-law spectrum f2 gw (/) = Q a f a ), the mathematical analysis of Sec. |V fj| goes 
through unchanged.^ The optimal way of combining the estimators 



1 



m " — — E (5-41) 



n v k=l 



for each detector pair is 



A := U ^r^ \ 11 , (5-42) 



whereQ 



K = °t = nij fe V 2 ■ (5.43) 
The inverse variance for the optimal estimator 



^optimal := fi\\ ii=cr - 2 (5.44) 



is given by 



36 Although cross-correlation signals ^'S and ( kl )S taken during the same time interval T are 
correlated, i.e., 

cov{ iij) S, (kl) S} := (^S (kl) S) - ( (i %)( (fc %) + , 

the variance-covariance matrix CVy-vj-n := cov{^S, is dominated by the diagonal terms 

%)(ii) = cov{MS,MS} =: in the large noise approximation. Thus, in practice, one can 

treat all of the measurements ( 12 ^S\, ^ 12 - ) 5 , 2, ••• , ^~ 1 '^S ni _ l p as effectively uncorrelated. This 
is because the detector noises (which are statistically independent of one another) are the only 
contributors to this approximation. 

37 For completeness, we note that the optimal way of combining correlated random variables 
xi,---,x m is given by jl := Ej=i Aj, where Aj = J2f=i and is tlie 

inverse of the variance-covariance matrix := cov{xi,Xj}. When the variance-covariance matrix 
is dominated by the diagonal terms Cu = cov{xi, Xi} =: af, the optimal combination of data 
reduces (approximately) to the uncorrelated result Xi = a~ 2 . As argued in the previous footnote, 
this is what happens for the cross-correlation signals Xi *-* W'S. 
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I I 



^optimal _ X X n *i ^ a 



(5.45) 



i=l j<i 



The squared signal-to-noise ratio for the optimal estimator is given by 



l i 



SNR^ timal = EE^ fo) SNR 2 



(5.46) 



i=l j<i 



These equations should be compared with Eqs. (|5.28|) and (|5.29|) in Sec. |VE . 

Using Eq. ( |5.45| ) for the inverse variance of the optimal estimator, we can derive an 
expression for the minimum value of fi after 4 months of observation (i.e., T tot = 10 7 sec), 
for a false alarm rate equal to 5% and a detection rate equal to 95%, where we optimally 
combine data from different detector pairs. In order to simplify the analysis, we will assume 



that Ti 



tot 



nT, where n :- 



Hi 



is the same for all pairs of detectors. 



For a single pair of detectors, the minimum value of S7 is given by 



( ii)Q 95%,5% = 2_L65 v77i;/) 



sec 



where 



IOtt 2 



df 



7 2 (|/|) 



(5.47) 



(5.48 



The above expression for ^a 2 follows from Eqs. ( |3.68| ) and ( j.73 ), together with the nor- 
malization condition fi = Q T for a stochastic background having a constant frequency 
spectrum Q gw (f) = Qq. The factor of 2 • 1.65 in Eq. ( |5.47| ) comes from the choice of the false 



alarm and detection rates. (See, e.g., Eq. (|4.37|) .) Table [Q| in Sec. [VI D| contains values of 



95%, 5% 




hf 00 for different detector pairs. 



For the optimal combination of data from multiple detector pairs, we have 



(n 



95%, 5% | 




I optimal 



2-1.65 

10 7 sec 



-2 



I I 



-1 



i=l j<i 



(5.49) 



This can also be written as 



i i 



K % ' 5% loptimal) =EE( (I,) C % ) • (5-50) 

i=l j<i 

Tables p^ - ^V| in Sec. |VI D| contain values of £7q 5% ' 5% | optimal for the optimal combination of 
data from multiple detector pairs, for all possible triples and quadruples of the five major 
interferometers. 
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D. 4-detector correlation 



Rather than combine data from multiple detector pairs as described in the previous 
section, one can directly correlate the outputs of 4 detectors (or, in general, 2N detectors) 
in a manner analogous to the single 2-detector correlation described in Sec. [Ill (Jl -p 8 ! In fact, 
it turns out that one can write down expressions for the optimally-filtered squared signal- 
to-noise ratio and the inverse square of f2o 5% ' 5% for the 4-detector correlation in terms of a 
simple sum of products of the corresponding quantities for the individual detector pairs. As 
noted above, the analysis given here can easily be generalized to the case of 2N detectors. 
At the end of the section, we summarize our results by writing down the key equations in 
the general 2A-detector form. 

To begin, we define the 4-detector correlation signal S to be the integrated product 



rl/Z rl/Z 

S := / dti--- dt 4 suti) • ■ • s 4 U 4 ) Q(ti, • ■ ■ , £4) , (5.51) 

J-T/2 J-T/2 

where 

Si(t) : = tn(t) + ni(t) (5.52) 

are the outputs of the detectors (i = 1, • • • , 4), and Q(ti, ■ ■ ■ , i 4 ) is an arbitrary filter function, 
which we will determine shortly. To save some writing in what follows, we will use the 
shorthand notation 





[T/2 


,T/2 


[ T/2 dH EE 


/ dh-- 


• / dti 


J-T/2 


/-T/2 


J-T/2 



(5.53) 

Throughout, we will assume that the gravitational strains hi(t) satisfy the statistical prop- 
erties listed in Sec. [EI B| . We will also assume that the noise rii(t) intrinsic to the detectors 
are (i) stationary, (ii) Gaussian, (iii) statistically independent of one another and of the 
gravitational strains, and (iv) much larger in magnitude than the gravitational strains. The 
goal is to determine the filter function Q(ti, ■ • ■ , i 4 ) that maximizes the signal-to-noise ratio 
SNR := n/a of S. 

Let us start by calculating the variance a 2 . Since we are assuming that the noise intrinsic 
to the detectors are statistically independent of one another and of the gravitational strains, 
and that they are much larger in magnitude than the gravitational strains, it follows that 

a 2 := (S 2 ) - (S) 2 w (S 2 ) (5.54) 
dH' (s 1 (t 1 ) ■ ■ ■ s 4 (t 4 )si(ti) • • ■ s 4 (t 4 )) Q(h, • ■ ■ , U)Q{t' l: •••,<) (5.55) 

dH' (n 1 (t 1 )---n 4 (U)n 1 (t' 1 )---n 4 (t' i )) Q{t±, ■ ■ ■ ,U)Q{t'i, • • • A) (5-56) 

dH' (ni(ti)ni(*i)) ■ ■ ■ (n 4 (t 4 )n 4 (t 4 )) Q(h, ■ ■ • , U)Q(t[, ■ ■ ■ , Q . (5.57) 
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[T/2 


-T/2 




'-T/2 



38 The correlation of 3 detectors (or, in general, 2iV + 1 detectors) yields a signal that has zero 
mean. 
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a 2 



(See Eqs. (|3.19| )- (l3~2~2| ).) Using the definition ( |3.23|) of the noise power spectra Pi(\f\), 
Eq. ( 5.57 ) can be rewritten as 

2) j_ J*f Pii\h\)---P*{\h\) 

r T/2 

/ dH ^hti...^hU Q( tu ---,U) 

J-T/2 
rT/2 

/ dH' e - i27r/l *i ■ • ■ e - i2nf ^ Q(t' ■ ■ ■ , Q , (5.58) 

J-T/2 



where 



00 ^00 
4 . 



—00 



d*f = / dh ■ ■ ■ / d/ 4 • (5-59) 



If we further define the (finite-time) Fourier transform 

rT/2 
-T/2 

which has as its inverse (for —T/2 < ti, ■ ■ ■ , £4 < T/2) 



/•T/2 

Wi, • • • , U) ■= / ^ e"**** • • • e- 12 ^ Q(h, ■ ■ ■ , t 4 ) , (5.60) 

J-T/2 



/oo _ 
d 4 /e&M ... e « 2 *Q (/l ... /4) ; (5 . 61) 
-00 



we obtain 

1 



° 2 w G) /I rf4/ Pi(|/i|) ■ ■ ■ P4(|/4|) ■ ■ ■ ' /4)|2 • (5 - 62) 

This can be written in an even more convenient form if we define an inner product 

POO 

{A,B):= / d i fP 1 (\f 1 \)---P±(\f4)A*(f 1 ,...,U)B(f 1 ,...,f A ), (5.63) 



where A(fi, ■ ■ ■ , f±) and B(fi, ■ ■ ■ , f^) are any two complex-valued functions of four variables. 
(See Eq. ( |3.69| ).) Note that the inner product is positive-definite since Pi(\f\) > 0. Using 
gTB5| ) it follows that 



o 2 »(1)\q,Q) • (5.64) 



4 

2 

To calculate the mean value //, we proceed in a similar manner. Since the noise intrinsic 
to the detectors are statistically independent of one another and of the gravitational strains, 
it follows that 

fi := (S) (5.65) 

rT/2 

= dH Mh) ■ ■ ■ Si {U)) Q{t ir --,U) (5.66) 

J-T/2 
T/2 

dH (h 1 {t 1 )---h±(U)) Q{h,---,t 4 ) . (5.67) 

-T/2 
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This can be expanded to 

/•T/2 

/i= / dH ({h 1 {t 1 )h 2 {t 2 )){h 3 {t 3 )h i {U)) + {13){2A) + {U){23)) Qfa,..., t 4 ) (5.68) 



-T/2 

by using the "factorization" property 

(£1X2X3X4) = (xix 2 )(x 3 x 4 ) + (xix 3 )(x 2 x 4 ) + (xix 4 )(x 2 x 3 ) (5.69) 

for Gaussian random variables xi,x 2 ,x 3 ,x 4 each having zero meanQ 

To express /x in terms of the inner product (|5.63|) , we first use (|5.61|) to rewrite 
Q(ti, ■ ■ ■ , t 4 ) in terms of its Fourier transform Q(fi, • ■ ■ , / 4 ): 



x 



T/2 

dH Uhtih^^ihsits^iU)) + (13)(24) + (14)(23) 

-T/2 v 

/DC _, 
tf fe i2nhh... e i*«uuQ( fir ..jj (5 70) 

-00 



Then by interchanging the order of integrations and rearranging terms in the integrand, we 
see that 



d 4 f / dH e i27Tfltl ■ ■ ■ e i27Tf * u 

-oo J-T/2 

x((/ il (t 1 )/ i2 (i 2 ))(/ i3 (t 3 )/i4(t4)) + (13)(24) + (14)(23)) Q(fi,---,U) , (5.71) 
or, equivalently, 

M = (AQ) , (5.72) 

where 

1 rT/2 

A(fu ■ ■ ■ , h) := -yrr^ / ci 4 t e"* 2 ^ 1 ■ ■ ■ e"* 2 ^ 

Pi(|/i|)---P4(|/ 4 |) 

x((/ii(fi)/i 2 (t 2 ))(/i 3 (t 3 )/i 4 (t4)) + (13)(24) + (14)(23)) . (5.73) 



To simplify this expression for A(f\, ■ ■ ■ , / 4 ), we expand the expectation values 
(hi(ti)hj(tj)) as 

/oo 

d/ c **/ft-*) H .. {f) > (5 74) 
-00 

where 

= l/r3 ^ w(l/l) 7 ^ (l/l) ' (5 - 75) 



39 In Eq. (13)(24) is used as a shorthand notation for {hi(tx)h 3 (t 3 )) fa^hifa)) , etc. 
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and Jij(f) denotes the overlap reduction function between the ij detector pair. (See 
Eqs. (ETSgp and (ggg) .) Substituting Eq. (|5Jl into ( jTTg) yields 



/4 



1 



Pi(\h\)---P4(\h\) 

r T/2 



/oo 
df H 12 (f)H u {f) 
-oo 



-T/2 

+13,24 + 14,23 , 



(5.76) 



where 13, 24 and 14, 23 denote the analogous terms with the appropriate interchange of 



detector indices 1, • • • , 4. Since 



T/2 
-T/2 



dt e i2wft 



(5.77) 



(see Eq. p.27|) , we can explicitly integrate over the time variables t\, ■ • ■ , t±: 

1 



/oo 
df H 12 (f)H u {f) 
-oo 



Pi(\h\) • ~ p*{\M) 

x5t(/i - f)S T (f 2 + f)5 T (f 3 - f')5 T (f 4 + /') 
+ 13,24+ 14,23 . 



(5.78) 



Also, since the noise power spectra Pi(\fi\) and filter function Q(fi, ■ ■ ■ , f±) are not expected 
to vary much over the support, 1/T, of 5r(fi — f) and 5t(/3 ~ f), we are justified in 
approximating these two finite-time delta functions by ordinary Dirac delta functions 5(j\ — 
f) and 5(fs — /'). This approximation allows us to eliminate the integrations over / and /', 
yielding 



A(ft 



Hn{fx) 



^(l/ii)^(|/ 2 |) 
+13,24+ 14, 23 . 



*r(/i + / 2 ) 



P*(\M)P*(\h\) 



*r(/ 3 + U) 



(5.79) 



Given the above expressions for <r 2 , /i, and A(fi, ■ ■ • , fi), it is now a simple matter to 
evaluate the squared signal-to-noise ratio of the 4-detector correlation S, and to determine 
the filter function Q(f\, ■ ■ ■ , fa) that maximizes this ratio. In terms of the inner product 
( |5.63| ), we have 



SNR 



a 2 ~ 



(q,q) 



(5.80) 



As we saw already in Sees. [Ill C| and |V A| , such a ratio of inner products is maximized by 
choosing 



Q(f 



fa)=X A(fa,---,fa) 



(5.81) 



where A is an arbitrary (real) overall normalization constant. For this choice of Q(fi, ■ ■ ■ , fa), 
the value of the squared signal-to-noise ratio is given by 
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SNR^ ptimal « 2 4 (A, A) . (5.82) 

To find an explicit expression for SNR„ ptimal , we substitute ( pT9|) into the RHS of ( ^82|) 
and expand the product of A(f±, ■ ■ ■ , / 4 ) with itself. This leads to nine different terms: three 
diagonal (i.e., squared) terms and six off-diagonal terms. The diagonal terms are given by 

d.agona. tenns = ? £ *, p //^ m f M + A) ^fa ^ + « 

+13,24+ 14,23. (5.83) 

A typical off-diagonal term is given by 

/oo ] 

x5 T (/ a + / 2 )*r(/ 3 + / 4 )*r(/i + / 3 )5t(/ 2 + / 4 ) • (5.84) 

Let us evaluate each of these terms separately. First, for the diagonal terms, by ap- 
proximating one of the two finite-time delta functions in each of the factors 5f.(/i + f 2 ) and 
<^t(/3 + fs) i n ( |5.83|) by ordinary Dirac delta functions, and by evaluating the others at 
fi + / 2 = and /3 + / 4 = 0, we eliminate two of the integrations, and introduce a factor of 
T 2 : 

# 2 2 (/) r M , H Uf) 



/oo Hit) f°° 



L (|/|)P 2 (|/|) ^-oc ' P 3 (|/1)P4(|/'|) 

+13,24 + 14,23 (5.85) 
/3ff 2 \ 4 2 y°° 7l 2 2 (l/l)^ w (l/l) y" „ ; 714(1/1)^(1/1) 

" Loo 1 / 6 Pi(|/|)P 2 (|/|) /-oo J /' 6 P 3 (|/1)P4(|/D 

+13,24 + 14,23 (5.86) 
= ^SNR 2 ^SNR 2 + ^SNR 2 ^SNR 2 + ^SNR 2 ( 23 >SNR 2 , (5.87) 

where ^^SNR 2 denotes the squared signal-to-noise ratios for the optimally-filtered cross- 
correlation signal for the ij detector pair, which we derived in Sec. |III C| . 

For the off-diagonal term (|5.84j ), we can again approximate 5t(/i + / 2 ) and Srifs + / 4 ) 
by ordinary Dirac delta functions to obtain 

poo poo j 

off-diagonal term « 2' / * / df 3 mfim \ fl \ mi ,\)P im ^ 



oo 



xiJ 12 (/ 1 )F 3 4(/3)^i3(/i)^ 2 4(-/i) 4(/i + /a) • (5-89) 

By further approximating one of the finite-time delta functions in 5^(f\ + / 3 ) by an ordinary 
Dirac delta function, and by evaluating the other at /i + / 3 = 0, we eliminate one more 
integration and introduce a factor of T: 

off diagonal term ~ 2 4 T f°° df H ^f) H ^f) H ^f) H ^f) f5 90l 

off dxagonal term ~ 2 T j_Jf ^(i/Dp^^fl/D^fl/i) " ( 5 - 90 ) 
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But since this term grow like T, while the diagonal terms grow like T 2 , we can (for large 
observation times) ignore the contribution of the off-diagonal terms to the optimal signal- 
to-noise ratio. The final result is thus[^] 

SNR^ ptimal « ^SNR 2 ( 34 )SNR 2 + ( 13 >SNR 2 ( 24 >SNR 2 + ( 14 >SNR 2 ^SNR 2 . (5.91) 

Note that the optimal signal-to-noise ratio for the 4-detector correlation is quadratic in the 
signal-to-noise ratios for the individual ij detector pairs. This is because the 4-detector 
correlation signal S, given by Eq. ( |5.51|) , is quartic in the detector outputs Sj(tj), while 
the cross-correlation signal values ™'S are quadratic in the detector outputs. (See, e.g., 

Eq. qpsp .) 

Given (|5.91 ), we can now ask the question: "What is the minimum value of Qq required to 
detect a stochastic gravity-wave signal 95% of the time, from data obtained via a 4-detector 
correlation experiment?" Since 



/r/j, 
dH s^h) ■ ■ ■ s 4 (U) Q{h,---,U) (5.92) 
-T/2 

d 4 f h(fi) ■ ■ ■ h(fa) Q*(fi,---, fa) , (5.93) 

-oo 

is, effectively, a sum (over fa,- ■ ■ , fa) of a large number of statistically independent random 
variables (products of the Fourier amplitudes Si(/i), which are correlated only when |/j — 
fj\ < 1/T), the central limit theorem guarantees that S will be well-approximated by a 
Gaussian random variable. Thus, the statistical analysis of Sec. [V| is valid for the 4-detector 
correlation as well. In particular, for a false alarm rate a = 0.05, and for a detection rate 
7 = 0.95, the minimum value of Qo for the 4-detector correlation is determined by setting 
the signal-to- noise ratio after n observation periods equal to 2 • 1.65 (see Eq. ( |4.36| )). For 
a stochastic background having a constant frequency spectrum f2 gw (/) = Q Q and for a 
total observation time T tot := nT, the squared signal-to- noise ratio for the optimally-filtered 
4-detector correlation signal S can be written as 

snr 2 « n* t 2 r df Mm r # iUlfl) 

optimal tot J ^ J /6 p i(|/|) p 2(|/|) J ^ J //6 p 3(|//|) p 4(|//|) 

+13,24 + 14,23. (5.94) 
Setting SNR optima i = 2-1.65 and rearranging terms yields 



^0 5% ' 5% Ltimal)" 4 (2 • L65) 2 - I -S I T t i t I df 



3#o 2 V^2 /■« 7^(1/1) 



°e timai ; v ' ; \10^J tot 7-oc J / 6 P 1 (|/|)P 2 (|/|) 

roo ^2 /I fl\\ 

X /-^ / '?™^ + 13 - 24 + 14 - 23 - (5 ' 95) 



40 This expression should be compared with 

SNR 2 ptimal = ( 12 )SNR 2 + ( 13 )SNR 2 + • • • + ( 34 )SNR 2 , 

which is the squared signal-to-noise ratio that we found in the previous section for the optimal 
combination of data from multiple detector pairs. (See Eq. ( 5.46| ).) 
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which can also be written in terms of the minimum values 0i)^ 5% ' 5% f or the individual ij 
detector pairs: 



(5.96) 



How does this minimum value of Qq compare with those found in previous sections for a 
single 2-detector correlation, and for the optimal combination of data from multiple detector 
pairs? First, from Eq. (|5.95 ) we see immediately that 



n 95%,5%| n-,-1/2 n7 x 

^0 (optimal ~ T tot • (5-97) 

This dependence on the total observation time T tot is the same as that for to')^ 5% > 5% f or a 
single detector pair, and for the optimal combination of data from multiple detector pairs: 

K % ' 5% loptimal) _2 = ( (12) ^ 5% ' 5% )" 2 + ( (13) ^ 5% ' 5 T 2 + " " " + ( (34) ^ 5% ' 5% )" 2 • ( 5 - 98 ) 

(See Eq. ( |5.50| ).) Thus, by correlating 2, 4 (or even 2N) detectors, one does not change the 
general dependence of the minimum value of Qq on the total observation time T tot . However, 
the numerical factors multiplying T tot differ from one another. In fact, it is fairly easy to 
show that ^o 5% ' 5% I optimal f° r the single 4-detector correlation is always greater than that 
for the optimal combination of data from multiple detector pairs. Thus, in theory, it is 
better to optimally combine data from multiple detector pairs than to optimally filter a 
single 4-detector correlation. Table [V| in Sec. [VI D| lists values of f^^ 5% ' 5% | optimal h\ m for the 
4-detector correlations taken from the five major interferometers: LIGO-WA, LIGO-LA, 
VIRGO, GEO-600, and TAMA-300. 

Finally, to conclude this section, we rewrite the key equations derived above for the 
general case of 2N detectors: 

(i) The 2iV-detector correlation signal S is defined by 

S := / d 2N t sxih) ■ ■ ■ s 2N (t 2N ) Q(t u • • • , t 2N ) . (5.99) 

J -T/2 

(ii) The "factorization" property for 2N Gaussian random variables each having zero mean 

(xix 2 ■ ■ -x 2 n) = (xix 2 )(x3X4) ■ ■ ■ (x 2 n-ix 2 n) + all possible permutations . (5.100) 



41 The number of terms on the RHS of the "factorization" equation is given by (2iV — 1) (2iV 
3) • • • 1. For N = 1, 2, 3, • • • , this corresponds to 1, 3, 15, • • • terms. 
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(iii) The inner product for which 



/ 1 \ 2N 

H=(A,Q) and ^=(2) ( 5 - 101 ) 

is given by 

/oo 
d 2N f Pi(|/i|) • ->P2„{\M) A*(f ± , • • • , f 2N )B(f 1 , ■■■J 2N ), (5.102) 
-00 

where A(fi, • • • , / 2 at) and • • • , /W) are any two complex-valued functions of 27V 

variables. 

(iv) The optimal filter function Q(t±, • • • , £ 2 tv) is given in the frequency domain by 

Q(fi,---j2 N ) = \A(f 1 ,-..J 2N ) , (5.103) 

where 

#2JV-l,2Jv(/2JV-l) r /, , x \ 

p n7 hp /if uMMr-i + /2vJ 

" 2JV-1 U J2N-1 1 J-T 2iV"U J2N \ ) 

+ all possible permutations . (5.104) 

(v) The squared signal-to-noise ratio for the optimally-filtered 2iV-detector correlation is 
given by 

SNR 2 ptimal « (12) SNR 2 (34) SNR 2 • • • ( 2iV_1 ' 2iV )SNR 2 + all possible permutations . 

(5.105) 

(vi) The minimum value of f2 required to detect a stochastic gravity-wave signal 95% of 
the time, with a false alarm rate equal to 5%, from data obtained via a 2iV-detector 
correlation experiment, is given by 

(^ 5% ' 5% | optimal ) _2 = (2- 1.65)^-^ [ (( 12 )^ 5% ' 5% )" 2 (( 34 )^ 5% ' 5% )" 2 

x • • • ((2iv-i,2iv) fi 95%,5%^ 2 + aU posgible permutations 1 ^ , (5.106) 

where (*?')^ 5% > 5% are the analogous quantities for the ij detector pairs. 
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E. Correlated detector noise 



We have shown how to carry out an experimental search for a stochastic background 
of gravitational radiation, by correlating the outputs of widely-separated detectors. Our 
analysis assumed that any correlation between the two outputs arises only from a stochastic 
gravity-wave background. In this section, we address the validity of this assumption, and 
look at possible sources of instrumental contamination that give rise to a correlated signal 
between the separated detectors. Any source of correlated environmental noise or interfer- 
ence in the separated detectors mimics the correlation arising from a stochastic background, 
so it is important to understand the order-of-magnitude effects of any potential sources of 
such correlated noise. 

This subject has already been considered in some detail, both in the the published paper 
of Christensen [[J and, in more detail, in Chapter 7 of his Ph.D. thesis |24j| . In the thesis 
work, the following sources of correlated detector noise are analyzed: 

(i) Seismic noise, whose effects on initial LIGO are minimal. 

(ii) Fluctuations in the residual gas (for two interferometers sharing a common vacuum 
system such as the LIGO-WA site). The effects on initial LIGO are minimal, and of 
course there is no effect for separated detectors that do not share a common vacuum 
envelope. 

(iii) Acoustic noise, whose effects on LIGO are minimal. 

(iv) Cosmic ray showers, whose effects on LIGO are minimal. 

(v) Magnetic field fluctuations, whose effects on initial LIGO might be significant. 

In this section, we derive a general formalism for calculating the effects of such correlated 
detector noise, and illustrate this for the case of correlated magnetic fields, which Christensen 
concluded would be the most significant source of correlated fluctuations at two widely 
separated sites. 

Although these ideas can be generalized to multiple site locations, for simplicity we 
present only the two-site case. Correlated noise in two detectors can be described by the 
cross-spectral function C(f) defined by 

1 f°° 

MtMO) = ~ / $ C(f) . (5.107) 

Z J —CO 

Because the LHS is real, the cross-spectrum satisfies C(f) = C*(— /). If this cross spectrum 
is non-vanishing, then it produces a correlation between the two detectors which mimics the 
effect of a stochastic gravity-wave background. 

It is straightforward to determine the point at which a non- vanishing cross spectrum C (/) 
will significantly interfere with a stochastic background search. Such interference will not 
take place if the correlation arising from the cross spectrum of detector noise is significantly 
less than that arising from the stochastic gravity-wave background: 

(ni(t)n 2 (t')> < (hi(t)h 2 (t')) for \t - t'\ < d/c , (5.108) 
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where d/c is the light travel time between the two sites. Making use of Eqs. ( |5.107| ) and 
( |5.74| ), we see that correlated detector noise will not impair a stochastic background search 
if 



l/r 3 n^d/i) 7(1/1) 



(5.109) 



over the range of frequencies included in the optimal filter Q(f). For example, for the 
initial LIGO detectors, this range of frequencies is from about 40 Hz to 300 Hz, and in a 
4 month search, the expected level of sensitivity is about fl h\ m ~ 10~ 6 , so the RHS is 
~ 3 x 10~ 49 h\ 0Q sec. Thus, in order that correlated sources of noise do not interfere with 
initial LIGO's 4- month stochastic background search, one must have 



\C(f)\ «3x 10" 49 h\ m sec for 40 Hz < / < 300 Hz 



(5.110) 



This limit can be stated in an interesting way. If we compare the allowable cross spectrum 
C(f) with the intrinsic noise power spectrum P(f) ~ 10 -45 sec in each detector, we see 
that they differ by about four orders of magnitude. Thus, in order that correlated noise 
sources do not interfere with a J^-month long stochastic background search for initial LIGO, 
the correlated sources of noise must not contribute more than 1% of the motion of the test 
masses, in the frequency range from 40 Hz to 300 Hz. 

Essentially the same limit on correlated noise can be written in another fashion. If we 
make use of Eq. (|3.52|) , the contribution of correlated noise to the expected mean of the 
signal can be written as 



/^"correlated noise 



1 -t 

2 



dfC(f)Q(f) 



(5.111) 



Requiring that this be smaller than the magnitude of the variance a 2 (see Eq. ( |3.68|) ) leads 
immediately to the condition that the correlated noise will not interfere with a stochastic 
background search over an observation time T if 



|c(/)l 2 « 



Tf 



(5.112) 



As before, the contribution of correlated noise to the motion of the interferometer must be 
smaller than the intrinsic detector noise motion by a factor of (T/) -1 / 2 . For an observation 
time T of 4 months {i.e., 10 7 sec) sec and frequencies / ~ 100 Hz, this gives the same 1% 
bound as before. 

One can also give a precise formula showing the effects of correlated noise sources on the 
expected value /i of the signal. Making use once again of the expected mean value of the 
signal (|3.52|) , we can express the ratio of the signal mean arising from the correlated noise 
to that arising from the stochastic background as: 



/^correlated noise 



/^■stochastic background 



IOtt 2 
3H! 



IZo dfC(f)Q(f) 



f-oodf\f\-* Ogwd/DTd/DW) 



(5.113) 



This formula allows us to precisely determine the effect of any correlated source of noise on 
the signal value starting from a model or measurement of the correlation spectrum C(f). 
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In particular, if we assume that the stochastic gravity-wave background has a constant 
frequency spectrum J7 gw (/) = Q Q , then 



/'■correlated noise 



/'■stochastic background 



^correlated noise limit ^ 114) 



■ 



where ^correlated noise limit is (by definition) the smallest value of f2 that can be observed 
before the effects of correlated instrument noise interfere with the measurement. 

The effects of a given source of correlated noise can be modeled more precisely. Here, 
we work through one example: the effects of correlated magnetic field fluctuations. Chris- 
tensen's work PJ24| concludes that these are the most likely environmental source of corre- 



lated noise between two sites. 

The LIGO interferometers use small magnets to steer and push the optical elements (i.e., 
mirrors and beam splitters). These magnets are an integral part of the Length Sensing and 
Control (LSC) system. Forces are applied to these magnets with electromagnetic coils, and 
they are part of the servo loop that uses modulation techniques to lock and monitor the 
path length difference between the test masses. External magnetic fields, present in the 
environment, exert forces on these magnets thus constituting one of the different sources of 
instrument noise. 

The external magnetic field is of particular concern, because it propagates at the speed 
of light, and therefore can give rise to correlations between sites on the time scale d/c. The 
magnetic field in the laboratory consists of two parts: a "local" part and a "global" part. 
The local part is the magnetic field arising from the instrumentation, wiring, and power 
lines within the lab; the global part comes from Schuman resonances of the earth and the 
ionosphere, and lightning strikes over the surface of the earth. The local part of the magnetic 
field will not correlate between widely-separated sites; the global part, however, is likely to 



be highly-correlated. Studies in the Caltech 40-meter prototype lab |25j have shown that 
the spectrum of ambient magnetic fields is of order 10 -7 Gauss/-\/Hz, but most of this field 
is local. Christensen reports in his thesis on a number of studies that show that the global 
magnetic fields are well-modeled above 20 Hz by a power-law spectrum: 

/ * \ -°- 88 

PM) = A {mr,J (5 - 115) 

with A ~ 1.2 X 10~ 17 Gauss 2 /Hz during magnetically noisy periods such as thunderstorms, 
and A ~ 1.8 x 10~ 19 Gauss 2 /Hz during magnetically quiet periods. Separate measurements 
have shown that these fields have a coherence of order r ~ 1/2 in the frequency range 
of interest, over widely-separated (almost antipodal!) points on the earth's surface. These 
highly-correlated global fields are the main concern here. 

In order to reduce the effects of external magnetic fields on the test masses and optics, 
the four magnets on each optic or test mass are arranged to approximately cancel both the 
dipole and quadrupole parts of the magnetic field. The magnitude of the resulting force on 
a test mass may be described by: 



A A 2 
eo + jei + -^e 2 + • • 



(5.116) 
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Here, \i is the magnetic dipole moment of one of the magnets, B is the magnitude of the 
ambient magnetic field, I is the length scale over which the magnetic field is varying in the 
vicinity of the test masses and optics, and A is the separation between the magnets on the 
test masses and optics. The quantities e ,ei,--- are the fractional difference between the 
dipole, quadrupole,- • ■ moments of the different magnets (which are not perfectly matched). 
For the initial LIGO detectors, the preliminary design values of these quantities are ap- 
proximately n = 0.11 Amp • meters 2 /c, eo = 0.05, £ = 6 cm, A = 15 cm. This leads to a 
force 

F = kB (5.117) 



on the test masses, with k ~ 0.1 dyne/Gauss |[25|| . This force accelerates the test mass, 
producing an equivalent strain 

^ ~ Ji^WL • < 5 - 118 > 

where M = 10 kg is the mass of the optic, and L = 4 km is the length of the arm. This 
gives rise to a cross spectrum 

C(/) = JpWIW PbU) • (5119) 



where r is the coherence. Evaluating the integral in fl5.113|) and ( |5.114j ) we find that the 
limits on Q Q are 

q f 10~ 7 during magnetically "noisy" times; (z,-\9n\ 

^'correlated noise limit i -, r -i n— Q 1 n u 05 yo.lZUJ 

L 1.5 x 10 during magnetically quiet times. 

Our conclusion is that for the initial LIGO design (where sensitivities are on the order of 
6 x 10~ 6 for 4 months of observation), magnetic field induced correlations are not a concern. 
However, for advanced LIGO (where sensitivities are on the order of 6 x 10~ n for 4 months 
of observation), the magnets must be eliminated from the design or they will they will 
significantly constrain the measurements of (or limits placed on) Qq. 



VI. NUMERICAL DATA 

This section consists of a series of graphs and tables containing numerical data for the 
five major interferometers: (i) Hanford, WA LIGO detector (LIGO-WA), (ii) Livingston, LA 
LIGO detector (LIGO-LA), (iii) VIRGO detector (VIRGO), (iv) GEO-600 detector (GEO- 
600), and (v) TAMA-300 detector (TAMA-300). This data was derived from published site 
location/orientation information and detector noise power spectra design goals using 



the stochastic background data analysis routines contained in GRASP [27|. (See Sec. |VT l 
for more information about the computer code that we wrote to perform these calculations.) 
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A. Noise power spectra 



Figures |i0| - [l"3| show the predicted noise power spectra for the initial and advanced LIGO 
detectors, and for the VIRGO, GEO-600, and TAMA-300 detectors. Figure displays all 
of the noise power spectra on a single graph. Figure |15| shows the predicted noise power 
spectra for the "enhanced" LIGO detectors, which track the projected performance of the 
LIGO detector design over the next decade. The data for the noise power spectra displayed 
in all of these figures were taken from the published design goals |26j . 



LIGO Noise Power Spectra 




FIG. 10. A log-log plot of the predicted noise power spectra for the initial and advanced LIGO 
detectors. 



VIRGO Noise Power Spectrum 
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FIG. 11. A log-log plot of the predicted noise power spectrum for the VIRGO detector. 
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GEO-600 Noise Power Spectrum 
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FIG. 12. A log-log plot of the predicted noise power spectrum for the GEO-600 detector. 



TAMA-300 Noise Power Spectrum 
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FIG. 13. A log-log plot of the predicted noise power spectrum for the TAMA-300 detector. 
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Noise Power Spectra 

(for the major interferometers) 
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FIG. 14. A log-log plot of the predicted noise power spectra for all the major interferometers. 
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FIG. 15. A log-log plot of the predicted noise power spectra for the "enhanced" LIGO detec- 
tors, showing the probable evolution of the detector design over the next decade. 
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B. Overlap reduction functions 
Figures |TE| - [2l}| show the overlap reduction functions 7(/) for different detector pairs. 



Overlap Reduction Function 

(LIGO-WA and the other major interferometers) 



TAMA-300 

- VIRGO 
GEO-600 
LIGO-LA 
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FIG. 16. The overlap reduction function 7(/) for the LIGO-WA detector and the other major 
interferometers. Note that the overlap reduction functions for the more distant detectors have their 
first zero at lower frequencies than those for the more nearby detectors. 



Overlap Reduction Function 

(LIGO-LA and the other major interferometers) 
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FIG. 17. The overlap reduction function j(f) for the LIGO-LA detector and the other major 
interferometers . 
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Overlap Reduction Function 

(VIRGO and the other major interferometers) 
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FIG. 18. The overlap reduction function 7(/) for the VIRGO detector and the other major 
interferometers. Note that the VIRGO and GEO-600 detectors are sensitive to almost orthogonal 
polarizations. 
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FIG. 19. The overlap reduction function 7(/) for the GEO-600 detector and the other major 
interferometers. Note that the VIRGO and GEO-600 detectors are sensitive to almost orthogonal 
polarizations. 
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Overlap Reduction Function 

(TAMA-300 and the other major interferometers) 
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FIG. 20. The overlap reduction function j(f) for the TAMA-300 detector and the other major 
interferometers . 
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C. Optimal filter functions 



Figures ^T]-§3 show the optimal filter functions Q(f) for different detector pairs, for a 
stochastic background having a constant frequency spectrum f2 gw (/) = Q . The optimal 
filter functions are normalized to have maximum magnitude of unity. 



Optimal Filter Function 

(initial LIGO-WA and initial LIGO-LA detectors) 
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FIG. 21. The optimal filter function Q(f) for the initial LIGO-WA and initial LIGO-LA 
detectors, for a stochastic background having a constant frequency spectrum Q gw (f) = Qq. The 
optimal filter function is normalized to have maximum magnitude of unity. 
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Optimal Filter Function 

(initial LIGO-WA and VIRGO detectors) 
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FIG. 22. The optimal filter function Q(f) for the initial LIGO-WA and VIRGO detectors, for 
a stochastic background having a constant frequency spectrum Q gw (f) = Qq. The optimal filter 
function is normalized to have maximum magnitude of unity. 



Optimal Filter Function 
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FIG. 23. The optimal filter function Q(f) for the initial LIGO-WA and GEO-600 detectors, 
for a stochastic background having a constant frequency spectrum Q gw (f ) = Qq. The optimal filter 
function is normalized to have maximum magnitude of unity. 
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Optimal Filter Function 



(initial LIGO-WA and TAMA-300 detectors) 
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FIG. 24. The optimal filter function Q(f) for the initial LIGO-WA and TAMA-300 detectors, 
for a stochastic background having a constant frequency spectrum figw(/) = ^0- The optimal filter 
function is normalized to have maximum magnitude of unity. 



Optimal Filter Function 

(initial LIGO-LA and VIRGO detectors) 



f' ! '', A 

■ 

il A 




WW 
^ Ml II y 





0.0 100.0 200.0 300.0 400.0 500.0 

f(Hz) 



FIG. 25. The optimal filter function Q(f) for the initial LIGO-LA and VIRGO detectors, for 
a stochastic background having a constant frequency spectrum f2 gw (/) = Qq. The optimal filter 
function is normalized to have maximum magnitude of unity. 
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Optimal Filter Function 

(initial LIGO-LA and GEO-600 detectors) 
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FIG. 26. The optimal filter function Q(f) for the initial LIGO-LA and GEO-600 detectors, for 
a stochastic background having a constant frequency spectrum fi gw (/) = Qq. The optimal filter 
function is normalized to have maximum magnitude of unity. 



Optimal Filter Function 

(initial LIGO-LA and TAMA-300 detectors) 
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FIG. 27. The optimal filter function Q(f) for the initial LIGO-LA and TAMA-300 detectors, 
for a stochastic background having a constant frequency spectrum f2 gw (/) = Qq. The optimal filter 
function is normalized to have maximum magnitude of unity. 
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Optimal Filter Function 

(VIRGO and GEO-600 detectors) 
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FIG. 28. The optimal filter function Q(f) for the VIRGO and GEO-600 detectors, for a 
stochastic background having a constant frequency spectrum fig^/) = Qq. The optimal filter 
function is normalized to have maximum magnitude of unity. 



Optimal Filter Function 
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FIG. 29. The optimal filter function Q(f) for the VIRGO and TAMA-300 detectors, for a 
stochastic background having a constant frequency spectrum fig^/) = Qq. The optimal filter 
function is normalized to have maximum magnitude of unity. 
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Optimal Filter Function 

(GEO-600 and TAMA-300 detectors) 
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FIG. 30. The optimal filter function Q(f) for the GEO-600 and TAMA-300 detectors, for 
a stochastic background having a constant frequency spectrum Q gw (/) = Qq. The optimal filter 
function is normalized to have maximum magnitude of unity. 
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D. Signal-to-noise ratios and sensitivities 



Table | contains the values of the theoretical signal-to- noise ratios after 4 months (i.e., 
10 7 sec) of observation, for the optimally-filtered cross- correlation signal between different 
detector pairs, for a stochastic background of gravitational radiation having a constant 
frequency spectrum il gw (/) = fi = 6 x 1CT 6 hi 00 . Table [TJ contains the minimum values of 
f^o ^100 f° r 4 months of observation, for a false alarm rate equal to 5%, and for a detection rate 
equal to 95%, for cross-correlation measurements between different detector pairs. Tables 
[TV| contain the minimum values Qq h\ QQ for 4 months of observation, for a false alarm rate 
equal to 5%, and for a detection rate equal to 95%, for the optimal combination of cross- 
correlation measurements between multiple detector pairs, taken from all possible triples 
and quadruples of the five major interferometers. Table [V] contains the minimum values of 
Qo h\ m for 4 months of observation, for a false alarm rate equal to 5%, and for a detection 
rate equal to 95%, for optimally-filtered 4-detector correlations. (Note: The calculation of 
the signal-to-noise ratios and minimum values of Qq assume that the magnitude of the noise 
intrinsic to the detectors is much larger than the stochastic gravity-wave background. This 
corresponds to Eq. ( |3.75| ) in the text.) 



TABLE I. Theoretical signal-to-noise ratios after 4 months of observation, for the opti- 
mally-filtered cross-correlation signal between different detector pairs, for a stochastic background 
of gravitational radiation having a constant frequency spectrum £\w(/) = = 6 x 10~ 6 h^ Q0 . 





LIGO-WA 


LIGO-LA 


VIRGO 


GEO-600 


TAMA-300 


LIGO-WA 




3.45 


1.74 


5.09 x 10- 1 


6.12 x 10" 2 


LIGO-LA 


3.45 




2.10 


7.66 x 10- 1 


9.16 x 10~ 2 


VIRGO 


1.74 


2.10 




1.56 


9.14 x 10~ 2 


GEO-600 


5.09 x 1CT 1 


7.66 x lO" 1 


1.56 




1.32 x 10~ 2 


TAMA-300 


6.12 x 1(T 2 


9.16 x 10 -2 


9.14 x 10~ 2 


1.32 x 10~ 2 





TABLE II. Minimum values of ^ioo f° r 4 months of observation, for a false alarm rate 
equal to 5%, and for a detection rate equal to 95%, for cross-correlation measurements between 
different detector pairs. 





LIGO-WA 


LIGO-LA 


VIRGO 


GEO-600 


TAMA-300 


LIGO-WA 




5.74 x 10~ 6 


1.14 x 10~ 5 


3.89 x 10~ 5 


3.24 x 10" 4 


LIGO-LA 


5.74 x 10~ 6 




9.45 x 10~ 6 


2.58 x 10~ 5 


2.16 x 10~ 4 


VIRGO 


1.14 x 10~ 5 


9.45 x 10~ 6 




1.27 x 10~ 5 


2.17 x 10~ 4 


GEO-600 


3.89 x 10" 5 


2.58 x 10~ 5 


1.27 x 10~ 5 




1.50 x 10~ 3 


TAMA-300 


3.24 x 10" 4 


2.16 x 10~ 4 


2.17 x 10~ 4 


1.50 x 10~ 3 
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TABLE III. Minimum values of f2o h\ QQ for 4 months of observation, for a false alarm rate 
equal to 5%, and for a detection rate equal to 95%, for the optimal combination of cross-correlation 
measurements between multiple detector pairs, taken from all possible triples of the five major 
interferometers, 

Detectors 95% ' 5% | 

.Detectors Ji I optimal n ioo 



LIGO-WA, LIGO-LA, VIRGO 


4.50 x 10" 


-6 


LIGO-WA, LIGO-LA, GEO-600 


5.55 x 10" 


-6 


LIGO-WA, LIGO-LA, TAMA-300 


5.74 x 10" 


-6 


LIGO-WA, VIRGO, GEO-600 


8.28 x 10" 


-6 


LIGO-WA, VIRGO, TAMA-300 


1.14 x 10" 


-5 


LIGO-WA, GEO-600, TAMA-300 


3.86 x 10" 


-5 


LIGO-LA, VIRGO, GEO-600 


7.27 x 10" 


-6 


LIGO-LA, VIRGO, TAMA-300 


9.43 x 10" 


-6 


LIGO-LA, GEO-600, TAMA-300 


2.57 x 10" 


-5 


VIRGO, GEO-600, TAMA-300 


1.27 x 10- 


-5 



TABLE IV. Minimum values of Qq /i^ 00 for 4 months of observation, for a false alarm rate 
equal to 5%, and for a detection rate equal to 95%, for the optimal combination of cross-correlation 
measurements between multiple detector pairs, taken from all possible quadruples of the five major 
interferometers . 

Detectors fl I optimal ^100 

LIGO-WA, LIGO-LA, VIRGO, GEO-600 4.17 x 10" 6 

LIGO-WA, LIGO-LA, VIRGO, TAMA-300 4.50 x 10~ 6 

LIGO-WA, LIGO-LA, GEO-600, TAMA-300 5.54 x 10~ 6 

LIGO-WA, VIRGO, GEO-600, TAMA-300 8.28 x 10~ 6 

LIGO-LA, VIRGO, GEO-600, TAMA-300 7.27 x 10~ 6 



TABLE V. Minimum values of f2o ^ioo f° r 4 months of observation, for a false alarm rate 
equal to 5%, and for a detection rate equal to 95%, for optimally-filtered 4-detector correlations. 

Dptprtnrs q 95% < 5% | 
uetectors u I optimal n ioo 

LIGO-WA, LIGO-LA, VIRGO, GEO-600 6.49 x 10" 6 

LIGO-WA, LIGO-LA, VIRGO, TAMA-300 2.51 x 10~ 5 

LIGO-WA, LIGO-LA, GEO-600, TAMA-300 5.44 x 10~ 5 

LIGO-WA, VIRGO, GEO-600, TAMA-300 4.68 x 10~ 5 

LIGO-LA, VIRGO, GEO-600, TAMA-300 3.84 x 10~ 5 
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VII. COMPUTER SIMULATION 



In Sees. pT|-|V|, we described the data analysis and optimal signal processing required 
for the detection of a stochastic background of gravitational radiation. This analysis was 
in terms of continuous time functions and their associated Fourier transforms. But, in 
reality, when one performs the actual data analysis, continuous time functions will be re- 
placed by discrete time-series and Fourier transforms by their discrete frequency counter- 
parts. The discrete data can then be processed by computer code that takes the appropriate 
FFTs, constructs the optimal filters, whitens and windows the data, etc. We have written 
a number of functions (in ANSI-C) to do precisely this. These functions constitute part 
of a general-purpose data analysis package for gravitational-wave detection, called GRASP 



(Gravitational Radiation Analysis & Simulation Package) [27]. In this section, we describe a 
computer simulation (made up of these functions) that mimics the generation and detection 
of a simulated stochastic gravity-wave signal in the presence of simulated detector noise. 
Documentation and further information about the code can be found in the GRASP user's 
manual. 



A. Purpose 

The main reason for writing the computer simulation was to verify many of the theoret- 
ical calculations that were derived in the previous sections. Specifically, we wanted to see if 
the theoretically predicted signal-to-noise ratio SNR, for a stochastic background of gravi- 
tational radiation having a constant frequency spectrum fi gw (/) = Q , would agree with an 
"experimentally" determined signal-to-noise ratio SNR produced by the simulation. If the 
theoretical and experimental signal-to-noise ratios agreed, we could be confident that the 
theoretical calculations were correct. If they did not agree, we would know that something — 
either a theoretical calculation or a technical issue related to the simulation itself — needed 
further investigation. 

The theoretical and experimental signal-to-noise ratios were said to be in agreement if 
the relative error defined by 



relative error : = 



SNR - SNR 



SNR 



(7.1) 



was less than (or approximately equal to) the inverse of the theoretical signal-to-noise ratio 
after n observation periods. This is the error that one would expect (approximately 68% 
of the time) if we approximate the sample variance a 2 by the true variance a 2 , and use the 
fact that the sample mean 

(7.2) 



42 The total observation time is Tt Q t := nT, where T is the duration of a single observation period. 
For our simulations, T = 3.2768 sec. 
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can itself be thought of as a random variable with mean \i and variance cr 2 /n. (Recall from 
Sec. fV] that Si, S 2 , ■ ■ • , S n are n statistically independent measurements of the optimally- 
filtered cross-correlation signal S, each associated with one of the n observation periods.) 
Thus, 



SNR := 



pL fJ, ± a/y/n fJL 1 



a 



a 



SNR± 



1 



n 



n 



(7.3) 



which implies 



relative error : = 



SNR - SNR 



SNR 



n SNR 



(7.4) 



This criterion was satisfied by our simulation runs for both the initial and advanced LIGO 
detector pairs. 



B. Flow chart 

A "flow chart" for the simulation is as follows: 

1. Input the parameters defining the simulation. This can be done either interactively or 
via #def ine statements at the beginning of the simulation program. These parameters 
are: 

(i) the site identification numbers for the two detectors. 

(ii) the number of time-series data points iV to be used when performing the data 
analysis (i.e., FFTs, cross-correlations, etc.). N should equal an integral power 
of 2. 

(iii) the sampling period At of the two detectors. (Note: T := NAt is the duration 
of a single observation period.) 

(iv) the constant Q that defines the stochastic background frequency spectrum: 

(v) the total number of runs n that make up the simulation. (Note: T tot := nT is the 
duration of the total observation period.) 

2. Using the site identification numbers, obtain site location and orientation information, 
and information about the noise power spectrum and whitening filter of each detector. 
(This information is contained in input data files.) 

3. Using the site location and orientation information, construct the overlap reduction 
function 7(/j) for the two detectors. (Here /j := i/(NAt) where % — 0, 1, • • • , N/2 — 1. 
By convention, we ignore the value of 7(/), or any other function of frequency, at or 
above the Nyquist critical frequency /Nyquist := 1/(2 At).) 
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4. Simulate the generation of a stochastic background of gravitational radiation having a 
constant frequency spectrum Q^f) = Q . This can be done in the frequency domain 
by using a random number generator to construct (complex-valued) Gaussian random 
variables hi(fi) and h 2 (fi) having zero mean and joint expectation values: 

(KUdhU))) = ^~2 T *V , (7.5) 

mifiMf))) = ^-2 T /r 3 , (7.6) 

(KifiMf))) = ^~2 T *U /r 3 l{fi) • (7.7) 

These are just the discrete frequency versions of Eq. fl3.56|) specialized to the case 
^gw(/) = £V Note that the Fourier amplitudes h\{fi) and h 2 (fi) fall-off like f~ 3 / 2 . 

5. Simulate the generation of the noise intrinsic to the detectors, using the information 
contained in the noise power spectrum data files. This can be done in the frequency 
domain by using a random number generator to construct (complex-valued) Gaussian 
random variables ni(fi) and n 2 {fi) having zero mean and joint expectation values: 

("K/OniC/i)) = \ T P x (/0 , (7.8) 

(nl(fi)n 2 (fi)) = \t Sa P 2 (fi) , (7.9) 
(hl(f i )h 2 (f 1 )) = 0. (7.10) 

These are just the discrete frequency versions of Eq. Q3.64j ). 

6. Construct the Fourier amplitudes 

Si(/i):=^i(/i)+ni(/<) , (7.H) 
52(/i):=M/i)+"2(/0 • (7-12) 



7. Whiten the data in the frequency domain by multiplying §i(fi) and s 2 (fi) by the 
frequency components W\{fi) and W 2 (fi) of the whitening filters of the two detectors: 

:=Si(/0 , (7.13) 

52(/i) := s 2 (/0 ^ 2 (/i) . (7.14) 

This multiplication in the frequency domain corresponds to the convolution of s\{ti) 
and W\{ti) (s 2 (ti) and W 2 {U)) in the time domain. (Note: The purpose of whitening 
the data is to reduce the dynamic range of the corresponding power spectra.) 

8. FFT di(fi) and o 2 (fi) into the time domain to obtain the corresponding time-series 
0\(ti) and o 2 {U). (Here U := iAt where i = 0, 1, • • • iV — 1.) 
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9. Repeat steps 4-8 to obtain another set of time-series data oi(tj) and 02(ti). Distinguish 
these two different sets of data with superscripts: ^Oi(ti), ^h 2 (ti), < - 2 - ) Oi(tj), ^o 2 (ti). 

10. Offset {1) Oi(ti) and {2) Oi(t;) ( (1) o 2 (^) and (2) o 2 (^)) by T/2, and combine them with one 
another and with data left over from the previous observation period to produce a 
continuous-in-time data set Oi(tj) (o 2 (ti)). Oi(£j) and o 2 (tj) represent the "raw" (i.e., 
whitened) signal+noise data streams output by the two detectors. 

(Note: Steps 4-10 make up the signal generation part of the simulation.) 

11. Test the input data oi(U) and o 2 {ti) to see if they have probability distributions con- 
sistent with that of a Gaussian random variable. If either set fails this test, reject 
them both, and repeat steps 4-10 to obtain new input data Oi(tj) and o 2 (ti). 

12. Window the data streams 0\(tj) and o 2 (tj) in the time domain, using a Hann window 
function to reduce side-lobe contamination of the corresponding power spectra. 

13. FFT the windowed data into the frequency domain to obtain the corresponding Fourier 
amplitudes 5i(/j) and 6 2 (/j). 

14. Unwhiten the data in the frequency domain by dividing di(fi) and d 2 (fi) by the fre- 
quency components Wi(fi) and W 2 (fi) of the whitening filters of the two detectors: 

h(fi) 
h(fi) 

15. Shift the input data streams 0\{ti) and o 2 {ti) forward in time by T/2, and repeats steps 
(12-14), obtaining another set of Fourier components si(/i) and s 2 (fi). Distinguish 
these two different sets of data with superscripts: ( - 1 - ) 5i(/i), ^s 2 (fi), ^Si(fi), ^§ 2 (fi)- 

16. Average (%(/,) and ^h(f t ) (^s 2 (f t ) and Mfi)) to produce ~ Sl (fi) (S 2 (/i)). §i(fi) 
and s 2 (fi) are the Fourier components of the unwhitened time-series data Si(tj) and 
S2(U)- (Note: The purpose of this averaging is to reduce the variance in the estimation 
of the spectra 5i(/j) and s 2 (fi)-) 

17. Construct the optimal filter function Q(fi) with the overall normalization constant A 
chosen so that /i — Qo T, using the noise power spectra specified by the input data 
files. 

18. From Si(/j), s 2 (fi), and Q(fi) calculate the optimally-filtered cross-correlation signal 
S corresponding to a single observation period T. 

(Note: Steps 11-18 make up the signal analysis part of the simulation.) 

19. Repeat steps 4-18 n times, generating a set of optimally-filtered cross-correlation signal 
values: Si, S 2 , ■ • • , S n . 



Ol(fi) 

Hfi) 
W 2 (fi) 



(7.15) 
(7.16) 
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20. From S\, S2, ■ ■ ■ , S n construct the sample mean 



and sample variance 



1 - 
/i := -22 Si 



(7.17) 



a 2 :-- 



n 



1 1=1 



(7.18) 



The sample (or "experimental") signal-to-noise ratio produced by the simulation is 
given by 



SNR := ^ . 
a 



(7.19) 



21. Calculate the theoretical signal-to-noise ratio, using a discrete frequency approxima- 
tion to the integral 



SNR = n 



10vr 2 



o / /Nyquist j f 7 2 (/) 
^ '0 7 / 6 Pi(/)P 2 (/) 



1/2 



(7.20) 



See Eq. (g?75|) . (Note: Since the dat a is discretely sampled, we should only integrate 
up to the Nyquist critical frequency /Nyquist := 1/(2 At).) 



22. From SNR and SNR, calculate the relative error 

SNR - SNR 

relative error : = 



SNR 



(7.21) 



As mentioned in Sec. |VII A| , this should be compared with the inverse of the theoretical 
signal-to-noise ratio after n observation periods l/(y/n SNR). 

(Note: Steps 20-22 make up the statistical analysis part of the simulation.) 

Note: In order to obtain signal-to-noise ratios on the order of 10 after n = 1600 runs, 
we needed to use rather large values of fi (e.g., 10~ 3 for the initial LIGO detectors, and 
10~ 8 for the advanced LIGO detectors). These large values meant that expression ( |7.20| ) 
for the theoretical signal-to-noise ratio had to be modified to properly take into account 
the contributions to the theoretical variance a 2 that are due to a large stochastic gravity- 
wave signal. (See Eqs. ( |5.3j ) and (p.4|).) Without these modifications, the theoretical and 
experimental signal-to- noise ratios would be more likely to disagree. Thus, instead of ( [T.20 ), 
we used a "mixed" expression for the theoretical signal-to-noise ratio: 



SNR 



T 



,/n r /Nyquist ir 7 2 (/) 

Jo aj f6 p l{f) p 2{f) 



J /Nyquist df 



pptu)p'iu) 



R(f) 



1/2 • 



(7.22) 



This involves Eqs. ( |5.3|) and ( |5.4|) for the variance a 2 , but the large noise expression (|3.73|) 
for the optimal filter function Q(f). 
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VIII. CONCLUSION 



In this paper, we derived the optimal signal processing strategy required for stochastic 
background searches. We discussed signal detection, parameter estimation, and sensitiv- 
ity levels from a frequentist point of view. We also discussed the complications that arise 
when one considers: (i) arbitrarily large stochastic backgrounds, (ii) non-stationary detector 
noise, (iii) multiple detector pairs, and (iv) correlated detector noise. We explained how we 
verified some of the theoretical calculations by writing a computer simulation that mimics 
the generation and detection of a simulated stochastic gravity-wave signal in the presence 
of simulated detector noise. And we noted that the "experimental" results and theoretical 
predictions agreed to within the expected error. These results suggest that both the the- 
oretical signal processing formulae and the implementation of these formulae in computer 
code are correct. But we should not stop here. For example, there are still a number of ways 
that we can improve the data analysis code before we use it to search for a real stochastic 
background in the outputs of real interferometers. To conclude this paper, we list some of 
the desired improvements below: 

(i) The first change that we would like to make is to calculate real-time noise power spec- 
tra for the detectors, and to use this calculated data (rather than the information contained 
in the input noise power spectrum data files) to construct the optimal filter function Q(fi). 
(See step 17 in the computer simulation described in Sec. |VIIB[ ) Since the real-time noise 



power spectra will change slightly from one measurement to the next, we could then apply 



the data analysis strategy discussed in Sec. VB for nonstationary detector noise 



(ii) In order to obtain accurate real-time noise power spectra for the two detectors, it will 
probably be necessary to use more sophisticated spectral estimation techniques. Currently, 
we use a Hann window to reduce side-lobe contamination, and we average two overlapped 
data sets to reduce the variance, when forming our estimates of Si(/i) and s 2 (fi). (See steps 
12, 15, and 16 of the computer simulation.) This procedure can be replaced by multitaper 
spectral estimation methods, which use a special set of window functions — called Slepian 
tapers — to form spectral estimates of time-series data.0 GRASP [27] contains a modified 



version of a public domain package by Lees and Park J5(J to perform the multitaper spectral 
estimation. In addition to providing better spectral estimates, multitaper methods also 
provide nice techniques for "spectral line" parameter estimation and removal. This feature 
will be extremely useful when analyzing data produced by a real detector. For example, 
one will be able to track contamination of a data set by the line harmonic at 300 Hz, and 
remove a pendulum resonance at say, 590 Hz. (See GRASP JSTJ for more information.) 



(iii) In addition to being able to identify and to remove "spectral lines" from a real data 
set, one would also like to be able to test the data to see if the distribution of sampled 
values is consistent with normal detector operation. For example, one might check the input 
data set to see if it has a probability distribution consistent with that of a Gaussian random 
variable. If the test reveals an exceptionally large number of "outlier" points, then that 



43 See the original paper by Thomson [gg] and the text by Percival and Walden [29 for more 
details. 
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particular data set can be rejected. (See step 11 of the computer simulation.) The GRASP 
data analysis package already contains a routine that performs this Gaussian test. But 
we would also like a more rigorously characterized test that compares the distribution of 
the current data with that during "normal" detector operation, which most likely is not 
Gaussian. 

(iv) Finally, even though it will still be a few years before we can analyze real data from 
any one of the major interferometers, real data from prototypes — like the Caltech 40-meter 
interferometer — can be used in computer simulations. For instance, rather than write a 
computer simulation (like the one we described in Sec. |V11| ) that mimics the generation and 
detection of a simulated stochastic gravity-wave signal in the presence of simulated detector 
noise, we can write a computer simulation that mimics the generation and detection of a 
simulated stochastic gravity-wave signal in the presence of real detector noise.0 The fact 
that the noise level of a prototype interferometer is much larger than that of a major inter- 
ferometer poses no problem; we can simply "dial-in" a larger stochastic background signal 
to be able to detect it in the same amount of observation time. Another nice feature of this 
fake stochastic background/ real detector noise simulation is that we can address all of the 
issues (i)-(iii) discussed above in a context where we can still compare "experimental" (i.e., 
simulation) performance against theoretical expectations. We must be totally convinced 
that the data analysis code is working as expected, before we can trust it when searching 
for a real stochastic background in the outputs of real interferometers. 

We are currently working on all of the above topics. We hope to complete and report on 
these projects sometime in the near future. 
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